t.>lo'.>ltx K run-{Way}.

2c .7?

L

E u...

 

 

 

 

 

 

 

 

 

am?“ Mr ,
s ..

u. «Fem ‘

34,33... : \

A».L.%l .h

«mm

.M ‘

C. . .
3 3...... ‘ it...
Q V 5-7 J I‘vw L | 0‘
‘ ‘5. ffimm
. was.

run. at:
figvflfi

I3

; i
..(r.

x. 15.3".”

1006

This is to certify that the
dissertation entitled

PARAMETER ESTIMATION IN
MULTI-AXIAL THERMAL DIFFUSIVITY EXPERIMENTS

presented by

Sean Edgar Davis

has been accepted towards fulfillment
of the requirements for the

Doctoral degree in Mechanical Engineering

 

”$44!”. Origlll/

Majod‘ o ssor’s Signature

1 J03. 1’0;

Date

MSU is an Affirmative Action/Equal Opportunity Institution

- _....-o-.-.-.-.-.--—.-.-._-._._ _._ V._ .

 

 

LIBRARY

Michigan State

University

 

 

 

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 c:/CIRC/DateDue.indd-p.15

PARAMETER ESTIMATION IN
MULTI-AXIAL THERMAL DIFFUSIVITY EXPERIMENTS

By

Sean Edgar Davis

A DISSERTATION

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

DOCTOR OF PHILOSOPHY
Department of Mechanical Engineering

2005

ABSTRACT

PARAMETER ESTIMATION IN
MULTI-AXIAL THERMAL DIFFUSIVITY EXPERIMENTS

By

Sean Edgar Davis
Thermomechanical analysis requires, in general, quantifying the thermal response of
materials and thus, the therrnophysical properties of thermal conductivity (or diffusivity)
and specific heat. Many materials yield an anisotropic thermal response, requiring
measurements of multiple components of the thermal conductivity tensor. The extended
flash method allows simultaneous measurement of multiple components of the thermal
diffusivity tensor along any axes of interest. The locations of the temperature sensors in
such an experiment have an affect on the ability to accurately estimate the desired
components of the diffusivity tensor. Prior studies have optimized the sensor geometry
for specific cases, but a general methodology for optimizing such an experiment while
allowing flexible positioning of the sensors during measurement has not been developed.
Here, D-optimization is applied to a simulated extended flash diffusivity experiment to
improve the accuracy of the extended flash method through optimization of the distance
between the sensors. The effects of experimental noise and heating power and duration
are also examined. Optimal sensor positioning was determined using two criteria:
maximizing the information obtained via all the estimated parameters, and maximizing
the amount of information obtained via the estimated diffusivities. Results indicate that
the optimal inter-sensor distance increases with an increasing ratio of in-plane to out-of-
plane diffusivity. The analytically determined optimal sensor positioning for an isotropic

material is validated via experimental measurements on A181 304 stainless steel, where it

is shown that the accuracy of the estimated parameters improves for data sampled at the

optimized locations.

When modeling the anisotropic thermal response of materials, the material may be
rotated such that the physical axes coincide with the principal axes of the thermal
diffusivity tensor, resulting in thermal orthotropy. During measurements of such a
tensor, however, the principal axes may be unknown, requiring a method to determine
principal values and the orientation of the principal directions while simultaneously
measuring the diffusivity. An analytical study was performed where the four non-zero
components of the diffusivity tensor a were estimated for a material possessing random
in-plane anisotropy on the order of certain manufactured or mechanically loaded
elastomers. Results indicate that a four-sensor array allows sufficient sampling of the
material response to permit estimation of a to within 1% of the reference values. When
orthotropy is assumed for a material exhibiting random in-plane anisotropy, the estimated
values of on, are resolved to within 0.4% of the reference values when the magnitude of
the anisotropy is on the order of that seen in common deformed or machined elastomers.
Should the orientation of the specimen be known prior to experimentation, it is shown
that sampling the temperature response in an orientation that is rotated from the principal
thermal axes may provide more information about the estimated parameters, thus

improving the accuracy of the experiment.

ACKNOWLEDGMENTS

Foremost, I would like to thank the members of my Ph.D. guidance committee, Drs. Neil
Wright, James Beck, Tom Pence, and Andre Lee, for their input and assistance during the
preparation of this dissertation. I would particularly like to thank Dr. Wright, my advisor,
for his encouragement and patience. I hope that his attention to detail, something we
both admire, is reflected in my preparation of this work. I am grateful for the support of
the National Science Foundation and the Whitaker Foundation, both of which generously
provided funding that made this research possible. I would like to thank Xin Huang for
her help in conducting experiments at Michigan State. I owe Dr. Edwin LeGall, a former
colleague in Dr. Wright’s lab, thanks for his encouragement in pursuing a Ph.D. I have
made many friends during my time here at Michigan State; these friends have provided
relief when difficulties were encountered in this research and their aid is gratefully
acknowledged. Finally, to my dear friends back in Baltimore who always made me feel

like I was home despite moving hundreds of miles away, I thank you for your support.

iv

TABLE OF CONTENTS

LIST OF TABLES .............................................................................. vi
LIST OF FIGURES ............................................................................ vii
NOMENCLATURE ............................................................................ x
CHAPTER 1. INTRODUCTION ............................................................... 1
CHAPTER 2. OPTIMAL POSITIONING OF SENSORS FOR THERMAL
DIFFUSIVITY MEASUREMENTS .......................................................... 7
2.0 Motivation .............................................................................. 7
2.1 Methods .............................................................................. 15
2.2 Results ............................................................................... 24
2.2.1 Measurement Duration ................................................ 24
2.2.2 Anisotropic Diffusivity ............................................... 28
2.2.3 Experimental Validation ............................................. 31
2.3 Discussion ........................................................................... 37
2.4 Closure ............................................................................... 45
CHAPTER 3. ON MEASURING THE THERMAL DIFFUSIVITY TENSOR
WHEN PRINCIPAL AXES ARE UNKNOWN A PRIOR! .............................. 46
3.0 Motivation ........................................................................... 46
3.1 Methods .............................................................................. 56
3.2 Results ................................................................................ 63
3.3 Discussion ........................................................................... 73
3.4 Closure ............................................................................... 79
APPENDICES .................................................................................... 80
Appendix A - Notes concerning IR thermography .................................. 81
Appendix B - Alternating Direction Implicit (ADI) program for modeling
orthotropic heat conduction in a slab specimen ............................. 88
Appendix C — Random in-plane anisotropy finite difference program ........... 128
REFERENCES .................................................................................. 152
v

LIST OF TABLES

Table 2.1 Estimated directional thermal diffusivity values and associated standard
deviations for a 1.50 mm stainless steel plate (average of four tests). 36

Table 2.2 Average deviation of the out-of—plane and in-plane diffusivities from the
reference value of CL = 3.95 x 10'6 m2/s for a 1.50 mm stainless steel plate. 36

Table 2.3 Values of A and A+ computed from measurements of the 1.50 mm steel

plate. 36
Table 2.4 Optimum sensor location based on A and A“ for various 7 in a phantom

elastomer with an out-of-plane diffusivity of 1.08x10'7 mz/s and tp = 1.0 s. 39
Table 2.5 Relation between t,,, optimum A and N, and variation in 0t. 40

Table 2.6 Relationship between noise level in measurements and the corresponding
accuracy in the thermal properties of interest for an isotropic specimen with
tp = 1.0 s. 42

Table 3.1 Representative values of B, on 1, (122, and an used in the analytical study. 61

Table 3.2 Error in estimating diffusivity when assuming orthotropy (i.e., neglecting
0112) for an anisotropic material with principal axes oriented by B with respect
to the measurement axes. The maximum shear diffusivity for this material is
about 18% of the median axial diffusivity value. 69

Table 3.3 Error in estimating diffusivity when assuming orthotropy for an anisotropic
material with principal axes oriented by B with respect to the measurement
axes. The maximum shear diffusivity for this material is about 35% of the
median axial diffusivity value for tests 18-22 and about 64% of the median
axial diffusivity value for tests 23-27. 69

vi

LIST OF FIGURES

Figure 2.1 Diagram of the top surface of the slab specimen illustrating the illumination
area and the in-plane geometry. The quarter plane modeled in the parameter
estimation scheme is bounded by the positive x- and y-axes, which themselves
are axes of symmetry. Typically, specimens are between 1.5 mm and 3.5 mm
thick. Three sensors are located on the bottom surface of the specimen along
the coordinate axes a distance d apart; the projections of these sensor positions
are shown on the modeled quadrant of the specimen.

Figure 2.2 Simulated temperature histories and inverse model reductions for an
isotropic specimen at sensor locations of d = 0, 7, and 20 mm: (a) tp = l s,
(b) tp = 10 s. O is defined as T-Ta, where T... is the ambient temperature.

Figure 2.3 Representative residuals, defined as the difference between the modeled
and measured temperature histories, for an isotropic specimen: (a) center
temperature sensor, (b) off-axis temperature sensor.

Figure 2.4 Panel (a) shows the in-plane temperature response and 11,- for an isotropic
specimen where d = 12 mm. The dots represent the in-plane temperature data,
the solid line the ax sensitivity coefficient, the dotted line the heat flux
sensitivity coefficient, the short dashed line the 0tz sensitivity coefficient, the
long dashed line the 0ty sensitivity coefficient, and the dot-dash line the
convection sensitivity coefficient. Panel (b) shows A (circles) and A+
(triangles) as a function of elapsed time for the same specimen and sensor
location.

Figure 2.5 In-plane temperature response (circles) and log A” (triangles) as functions
of time for an isotropic specimen with d = 12 mm.

Figure 2.6 Variation of A and A“ with respect to d for an isotropic specimen with
t,, = 1.0 5. Calculated values of A are represented by circles, while calculated
values of A” are represented by triangles.

Figure 2.7 Location of the maximum A (circles) and A+ (triangles) as a function
of 7.

Figure 2.8 Variation of A and A+ with respect to d for an isotropic steel specimen

with t,, = 1.0 5. Calculated values of A are represented by circles, while
calculated values of A+ are represented by triangles.

vii

16

25

26

29

3O

32

33

35

Figure 3.1 Isotherms on the bottom surface of a specimen at t = 30 5 following 1 s
heating. (a) transversely isotropic: an = an = 0.16 mmzls, on; = 0, (b)
anisotropic: om = an = 0.16 mmz/s, an = 0.0282 mmzls, (c) orthotropic:
a1. = 0.188 mmzls, an = 0.132 mmZ/s, on; = 0. Note that Figure 3.1c is
congruent to Figure 3.1b if rotated 45° counterclockwise. 49

Figure 3.2 Schema of the top surface of the specimen illustrating the illumination
area. Top and bottom faces are subject to specified constant and uniform
convection coefficients. Heating on the top face is of known power and
duration; lateral edges are adiabatic. The temperature sensors are located on
the bottom surface of the specimen at the comers of a 9 mm square that is
aligned with the coordinate axes and has one sensor at the center of the bottom
face; the projection of these sensor positions are shown on the top surface of
the specimen. The principal axes of a, x’ and y’, are oriented at an angle B
from the experimental axes; B is assumed unknown at the start of parameter
estimation. 58

Figure 3.3 Simulated temperature excess, (-3 = T-Tm, and inverse solutions for (a)
orthotropic test (case 1) and (b) maximum shear test (case 5). Note that for
case 5, the x- and y-sensor temperature histories are identical as a result of an
and (1.22 being equal. 64

Figure 3.4 Scaled sensitivity coefficients for case 1: (a) y- sensor, (b) comer sensor,
(0) center sensor, ((1) x-sensor. The solid line represents F for (133, the dotted
line F for on 1, the dashed line F for an, and the dot-dash line F for an. Time
is not scaled, as the choice of appropriate scaling parameters is not immediately
obvious. For example, if the specimen thickness and corresponding axial
diffusivity cm are used as scaling parameters, than 1 second would correspond
to a Fourier number of 0.0104; however, if either of the principal in-plane
diffusivities are used with the specimen width as scaling parameters, 1 second
would correspond to a Fourier number of either 8.23 x 10'5 or 1.18 x 104. 65

Figure 3.5 Scaled sensitivity coefficients for case 5: (a) y- sensor, (b) comer sensor,
(0) center sensor, ((1) x-sensor. The solid line represents F for cm, the dotted
line F for on 1, the dashed line F for an, and the dot-dash line F for (x12. 66

Figure 3.6 Input data and model fits at the x-, y-, and comer sensors for case 5
(B = 45°) where orthotropic behavior is assumed during parameter estimation.
The solid lines represent the model fits for the simulated temperature history
at the indicated sensor. The poor fit between the data and the model
temperature histories at the corner sensor is due to the assumption of no
shearing diffusivity. The temperature responses at the x- and y-sensor
locations overlap since the x- and y-diffusivities are equal. 71

viii

Figure 3.7 Normalized A as a function of measurement duration for cases 1 through 5.
Empty circles represent data taken for B = 0°, triangles B = 15°, squares B =
22.5°, diamonds B = 30°, and filled circles B = 45°. 72

Figure 3.8 Variation of A for an (I having principal values of [1.10, 1.88, 1.32] x 10'7
mzls. The amount of information obtained from a measurement increases as
the specimen orientation shifts away from the principal axes of a. Note that for
B = 0° and 25°, A is six and two orders of magnitude smaller, respectively, than
A at B = 45°. 78

ix

FOi
hb

h:

Pi

qo

NOMENCLATURE

specific heat capacity (J/kg K)

intersensor distance (mm)

Fourier number in the i-direction (= (It/[2) (ND)
bottom surface convection coefficient (W/m2 K)

top surface convection coefficient (W/m2 K)

identity matrix (ND)

thermal conductivity tensor (W/m K)

component of the thermal conductivity tensor in the ij-direction (W/m K)
halfwidth of heated area (m)

halfwidth of specimen (m)

number of sensors

number of data points

estimated parameter (varies)

heat flux (W/mz)

temperature (°C)

glass transition temperature (°C)

temperature measured at time j (°C)

ambient temperature (°C)

the maximum rear surface temperature excursion (°C)
time (s)

heating duration (3)

experiment duration (5)
thermal diffusivity tensor (m2/s)
thermal diffusivity component in the i-direction (mZ/s)

angle between the physical (experimental) and principal thermal axes of 0t
(degrees)

scaled sensitivity coefficient at sensor position i with respect to estimated
parameter j (K)

anisotropy ratio (ND)

specimen thickness (m)

Fisher Information Matrix (ND)

density (kg/m3)

standard deviation of thermal diffusivity (mz/s)
temperature excursion, T-T... (K)

chi-squared error (K2)

determinant of the Fisher Information Matrix (ND)

ratio of the determinant of the Fisher Information Matrix for all estimated
parameters to the determinant of the Fisher Information Matrix for the nuisance
parameters (ND)

del operator (m")

xi

CHAPTER 1

INTRODUCTION

As computational technology advances in fields such as engineering and medicine,
quantifying material thermophysical properties becomes increasingly important. The
response of many materials may be assumed to be isotropic; however, others such as
tissues, composites and polymers may exhibit anisotropic responses. This anisotropy
could be due to the material structure (e.g., crystalline or composite), manufacturing
processes (Doss and Wright, 2000) or deformation (LeGall, 2001). Doss and Wright
(2000), among others, have demonstrated that the thermal diffusivity in seemingly

isotropic polymers may vary by 30% or more.

While temperature is a scalar, thermal conductivity is, in general, a tensor. For such a
material, the three-dimensional heat diffusion equation may be written as

pc%f—=V(kVT) (1.1)

where p is the density of the material in question, c the specific heat, T temperature, I
time, and k the second-order thermal conductivity tensor. The thermal conductivity
tensor has nine components

k1 1 kn k13

k: k21 k22 k23 (1.2)
k31 [<32 1‘33

however, since k is a symmetric tensor (which follows from classical thermodynamics
and Onsager’s principle of microscopic reversibility (Onsager, 1931)), only six of the
components are independent (e.g., k23 = k32) (Powers, 2004). Commonly, k is assumed to
be isotropic, with k = k-I. In general, k may be a function of T; however, here only small
temperature changes are considered such that k is a function of direction only. For many
circumstances, anisotropy is considered such that the principal axes of k are coincident
with the physical axes of the problem; then, only k; 1, kn and k33 have nonzero values.
The principal axes of the thermal response are unknown for materials of unknown or
complex composition, requiring the off-axis terms to be considered. The thermal
diffusivity a, where a = (l/pc)-k, is often used instead of k in transient problems, and is

often easier to measure.

Possible reasons for anisotropy in polymers include reorientation of the molecular chains
or the scission and reforming of inter-chain bonds during loading of elastomeric materials
(Wineman and Min, 2003). Manufacturing processes of many plastics and elastomers
create anisotropy due to molecular reorientation during molding, calendaring, pressing,
and similar processes. Kilian and Pietralla (1978) illustrated thermal anisotropy in
uniaxially drawn polyethylenes, while Broennan et al. (1999) measured anisotropic
conduction in stretched silicone rubber. LeGall (2001) has shown that the thermal
diffusivity of certain elastomers change during finite biaxial deformation. The
orientation of the principal axes of thermal diffusivity may be changing during the
deformation of a material, especially during general deformations, thus it cannot be

assumed, in general, that the principal axes of conduction coincide with the principal axes

of deformation. The challenge is to develop a procedure that provides optimal
assessment of the magnitude of thermophysical properties of materials, as well as the

direction of the principal axes of these properties.

Methods to measure thermophysical properties must be robust enough to deliver accurate
and precise results while providing indication of the accuracy of the model on which they
are based. The estimation of thermophysical properties for a material involves two
processes: a physical experiment in which a material specimen is heated and temperature
history data is recorded, and a parameter estimation routine where the data are matched to
an analytical model of the physical experiment, allowing the determination of the
specimen material properties, often through an optimization scheme. Traditional
experimental methods of thermophysical property measurement can be divided into those
for thermal conductivity, usually steady state (e.g., the guarded hot plate; Bradley, 1943),
and those for thermal diffusivity, usually transient (e.g., the flash method; Parker et al.,
1961). Many methods require excellent contact with the material specimen and, thus, are
unsuitable for tests during which the specimens may be deformed. Contact with physical
temperature sensors such as thermocouples disturbs the temperature profile of the
specimen, further complicating the problem (Keltner and Beck, 1983). The self-heated
thermistor, originally proposed by Chato (1968), has been primarily used to measure
thermal conductivity. Valvano and others have extended its use for thermal diffusivity
measurements as well as perfusion in biological materials (Valvano et al., 1985; Valvano
and Nho, 1991; Yuan et al., 1993), but the derived properties are volume averaged and

thus do not allow measurement of thermal properties along axes of interest. Conversely,

non-contact thermo-optical methods are often used to measure directional components of
the thermal diffusivity tensor. These methods include the flash, thermal wave, and
transient gradient methods. While the original flash method, developed by Parker et al.
(1961) and standardized by ASTM (E 1461-92), only measures the component of thermal
diffusivity perpendicular to the plane of heating, extensions of the flash method allow
simultaneous measurements of diffusivity along any axes of interest (e.g., Maillet er al.,
1990; Graham et al., 1999) and may even permit the measurement of off-axis diffusivity
terms as well. A fully populated a tensor may be thus resolved with a single test,

assuming a properly designed experiment.

The parameters that are calculated in a parameter estimation routine may be broken down
into two categories: the parameters of interest, such as thermal diffusivity, conductivity,
or specific heat, and nuisance parameters, which are parameters that are not of interest to
the experimenter yet must be estimated in order to arrive at accurate estimates for the
parameters of interest (Beck and Arnold, 1977). Typically, nuisance parameters involve
boundary conditions and may include an imposed heat flux, convection coefficients, or
emissivities of high-temperature bodies. Ideally, the experimental design should render
solutions for the parameters of interest relatively insensitive to errors in either the
nuisance parameters or the geometric parameters of the experiment, such as sensor
location. Only reasonably accurate estimates of the nuisance parameters may be needed
in order to determine the parameters of interest to a high degree of confidence, thus, less
computational effort may be expended in finding the nuisance parameters (Beck and

Arnold, 1977). Location of the temperature sensors, however, often plays a significant

role in the accuracy of the experiment (for an illustration of the relative effect of sensor
positioning on the accuracy of thermal property estimation, see Mzali et al., 2004). The
selection of temperature sensor location, specimen size, and heating geometry and
duration a priori requires knowledge of the thermophysical properties to be measured.
Of these variables, however, only temperature sensor location may be optimized along
with other parameters post-experiment by measuring the temperature field history. The
use of geometrically fixed sensors prohibits redesigning the physical experiment “on the
fly”, so each time the optimization criteria indicate that the sensor locations should be
moved, the experiment must be performed anew with the appropriate adjustments made
to the sensor locations. Optical methods of temperature measurement such as focal plane
array IR cameras can overcome this limitation by allowing flexibility in the choice of
sensor position. The temperature sensitivity of available cameras is on the order of that
of typical commercial thermocouples, and sampling rates now typically exceed 60 Hz for
320 x 256-pixel detectors. Most importantly, IR cameras allow full-field measurements
on the specimen surface, although such field measurements yield large data files. Adding
this level of optimization to the thermophysical property estimation routine increases the
complexity of parameter estimation, but careful design of the experimental boundary
conditions and the data reduction model can yield methods whose reduction time is
comparable to previous non-optimized models for thermal property estimation. The use
of such full-field measurements allows the experimenter to redesign the experiment to
make best use of the information available to capture unknown components or
orientations of the thermal diffusivity tensor 0. Presented here is (a) an optimization of

sensor position and heating magnitude and duration in an analytical extended flash

diffusivity experiment performed upon a material exhibiting transverse isotropy, (b)
experimental verification of the analytical procedure using a steel plate, (c) an
investigation of the ability to measure the four nonzero components of the a tensor in a
specimen exhibiting monoclinic (in-plane) anisotropy where the orientation of the
principal axes of a is assumed unknown prior to testing, and (d) quantification of the
error inherent in the assumption of orthotropic behavior for a specimen that actually
behaves in a monoclinically anisotropic fashion, i.e., the orientation of the principal
thermal axes in one plane (where the behavior may be characterized as orthotropic) is not

aligned with the physical axes of the experiment.

CHAPTER 2

OPTIMAL POSITIONING OF SENSORS FOR

THERMAL DIFFUSIVITY MEASUREMENTS

2.0 Motivation

Advances in thermomechnical analysis of polymers, composite materials, and biological
tissues are leading to improved designs in manufacturing and medicine.
Thermomechanical analysis requires, in general, quantifying the thermal response of
materials and thus, the thermophysical properties of thermal conductivity (or diffusivity)
and specific heat. Many materials yield an anisotropic thermal response, requiring
measurements of multiple components of the thermal conductivity tensor. The
anisotropy may result from the microstructure of the material (e. g., crystallinity, grain
boundaries, molecular alignment in polymers) or the macrostructure (e.g., composite
fibers, wood grains). Calendaring of PVC, for example, can produce differences of 30%
in the components of thermal diffusivity of an apparently isotropic material (Doss and
Wright, 2000). Measurement of the thermal conductivity tensor may also provide a
nondestructive test method for the mechanical response of a material if the
thermophysical properties can be correlated to deformation (Choy er al., 1978; LeGall,

2000) or stress (Venerus er al., 1999).

Several thermo-optical methods are available to measure the components of the thermal
diffusivity tensor including thermal wave (Mandelis, 1985), transient grating (Venerus et
al., 1999b), and extended flash (Maillet et al., 1990; Doss and Wright, 2000). Each of
these methods involves measuring the temperature of a specimen in response to localized
heating. Design parameters include the specimen size and shape, the heating area,
duration, and magnitude, and the position and number of the temperature sensors, each of
which may be a function of the thermophysical properties themselves. Measuring
thermophysical properties usually involves estimating two types of parameters: the
parameters of interest, such as thermal diffusivity, conductivity, or specific heat, and
nuisance parameters, which are parameters that are not of interest to the experimenter but
must be estimated in order to arrive at accurate values for the parameters of interest
(Beck and Arnold, 1977). Examples of nuisance parameters may include the heat input
or convection coefficients. Ideally, experimental design should minimize the influence of
errors in either the nuisance parameters or the design parameters of the experiment. The
nuisance parameters may be estimated at a reduced accuracy with little effect on the
parameters of interest; thus, less computational effort may be expended in finding the
nuisance parameters (Beck and Arnold, 1977). Location of the temperature sensors,
however, often plays a significant role in the accuracy of the measurements (Emery and
Fadale, 1997). Calculating the optimal temperature sensor location, specimen size, or
heating geometry requires knowledge of the thermophysical properties to be measured.
Of these, only the location of the temperature sensors may be optimized along with the
estimated parameters by measuring the temperature field history using, for example, a

focal plane array IR camera. The use of such a thermo-optical device allows full-field

temperature histories to be recorded, and thus, choice of the sensor locations after the

experiment.

Recent studies of estimating thermal properties in multiple directions have not fully
considered concurrent optimization of the sensor position along with the estimated
thermal parameters. Maillet et al. ( 1990) extended the flash method (Parker et al., 1961)
to include simultaneous estimates of the radial and axial thermal diffusivities in a
cylindrical specimen subject to an infinitesimal heat impulse. They used a geometric
method to optimize experimental parameters such as heating radius and eccentric sensor
position. Because they used point temperature measurements, they could optimize
position only prior to measurements. Doss and Wright (2000) extended the flash method
to measure the thermal diffusivity of an anisotropic PVC sheet in three orthogonal
directions using an instantaneous heat impulse. LeGall (2001) substituted a finite-
duration infrared laser pulse for the instantaneous flash and measured the thermal
diffusivity in three orthogonal directions of several elastomers subjected to biaxial finite
deformation. Neither Doss and Wright nor LeGall, however, focused on optimizing the
experimental parameters or sensor positioning. Instead, the design of their fixed-position
thermocouple probe was based on the work done by Maillet et al. and not reoptimized for
the geometry or material properties under investigation. Petit and coworkers have
focused on replacing detailed numerical models of multidimensional finite flash
experiments with lower-order reduced models (Petit and Hatchette, 1998; Le Niliot et al.,
2000; Videcoq and Petit, 2001). These studies were primarily analytical and dealt with

finding unknown boundary or interface conditions, as opposed to determining

thermophysical properties. Graham and coworkers made 2-D thermal diffusivity
measurements in orthotropic composite materials using an IR camera to record
temperature histories of specimens at multiple locations (Graham et al., 1999; Graham et
al., 1999b). Their data reduction was based on a Green’s function formulation that
assumes three mutually orthogonal thermal diffusivities and provided simple temperature
ratio relations for estimating the diffusivities in two directions simultaneously. While the
position of the eccentric temperature sensors were estimated by Graham et al. in
computational experiments, the sensor positions themselves were not optimized, and they
examined only instantaneous heating of the specimen, which was shown unsuitable for
materials of low thermal diffusivity, such as polymers and biological tissue, due to large
temperature excursions on the top face and low temperature signals on the bottom face
(LeGall and Wright, 2000). Graham et al. also proposed a method for making
multidimensional thermal diffusivity measurements in materials with random in-plane
orthotropy, analogous to the strain rosette, but reported no such measurements. Emery
and Fadale (1997) demonstrated that imprecision in the location of temperature sensors
can have a significant effect on the estimated parameters. Sensitivity coefficients (which
describe the change in temperature response of a specimen with respect to a given
parameter) were computed for an analytical model of a slab experiencing convective
boundary conditions on both surfaces. Using the sensitivity coefficients, the determinant
of the Fisher Information Matrix, which serves as a measure of the confidence in the
estimated parameters, was calculated. Position uncertainty was shown to have the
greatest effect on the value of the determinant, and thus, the ability to estimate the

thermal properties when sensors were located close to points of high heat flux.

10

Appropriate design of an experiment can reduce the effects of uncertainty in sensor
positioning by indicating where points of maximum sensitivity lie via calculation of the
sensitivity coefficients and covariance of the estimated parameters over a range of

possible experimental designs (Fadale et al., 1995; Emery et al., 2000).

Several studies have examined optimal designs for measuring thermophysical properties
by analytical means, without experimental comparison. Sawaf and Ozisik (1995)
examined measuring orthotropic conductivity in a rectangular parallelepiped by assuming
that three sides of the sample are subject to a uniform, constant heat flux while the other
three sides are assumed adiabatic. They simulated an experimental temperature response
by computing the solution to the diffusion equation using hypothetical properties while
adding random noise (on the order of 0 to 1 K, varied between simulations) to the
calculated temperatures. They showed reliable convergence to the thermal conductivity
used to generate the simulated experimental response using a Levenberg-Marquardt
method with an error of less than or equal to 3.2% in each direction. Other recent studies
have examined the optimal physical dimensions, heating area, or sensor positions for
reliable property estimation. Mejias et al. (1999; 2000; 2003) studied the effects of
sensor positioning, heating time, experiment duration, and the use of insulated or constant
temperature boundary conditions on the estimation of three thermal conductivity
components for a parallelpiped specimen experiencing a constant flux over three
surfaces. Dowding and Blackwell (1999) studied the effects of sensor positioning,
heating time, heating magnitude, and cooling time on the ability to estimate thermal

conductivity and volumetric heat capacity on a l-D model. Aviles-Ramos (2002) studied

11

the sensitivity coefficients developed for the thermal conductivities and heat capacities
obtained from a specimen having an orthotropic layer overlaying an isotropic layer.
These studies, primarily analytical in nature, often incorporate idealized boundary
conditions that may be difficult to achieve physically. For example, Mejias er al. (1999)
specified the three heated faces to experience distinct uniform heat fluxes that are
functions of time only. Even if the output of the heating device is well characterized, it
may be difficult to ascertain precisely the energy that is absorbed by a specimen, or even
where that absorption occurs in some materials (see McMasters et al., 1999; LeGall,

2001), and heating multiple faces may lead to additional nuisance parameters.

Several studies have investigated the experimental validity of optimal design. The
graphical analysis by Maillet et al. (1990) was one of the first attempts to optimize both
sensor positioning and heating geometry for multiaxial thermal property measurements.
Taktak (1992) investigated sensor position and variation of the heating time and
experimental duration on the optimality of measurements of thermal conductivity in
composites based on the guarded hot plate method. Taktak used D-optimality wherein
the determinant of a matrix of nondimensional sensitivity coefficients of the parameters
of interest (i.e., effective conductivity and volumetric heat capacity for Taktak et al.,
1993) served as the basis for comparison. Maximizing this determinant results in the
minimization of the variation in the parameters of interest. Theoretical results based on
l-D geometries were verified using a specimen containing embedded thermocouples.
McMasters (1997) considered several different models for the 1-D heating of a specimen

that involved departures from the classical flash experiment (i.e., internal radiation,

12

penetration of the flux boundary condition into the specimen, variable conductivity, and a
two-sided flash accounting for reflection to the bottom of the specimen surface).
Residuals, defined as the difference between the modeled and measured temperatures,
were examined to demonstrate the accuracy of the analytical model, while sensitivity
coefficients were plotted to illustrate the relative independence of the estimated
parameters. McMasters focused on optimizing the data analysis rather than the
experimental design. Mandelis and coworkers have developed a thermal wave method of
thermal diffusivity measurement that uses the thermal response of the heated surface to
determine the thermal characteristics of the specimen. The theory for the thermal wave
method has been well described by Mandelis (1985), but the technique may be better
suited for detection of cracks or inhomogeneities in materials, rather than measuring
thermal diffusivity (Mandelis et al., 1996). Mzali et al. (2004) analyzed 2-D flash
diffusivity measurements in which an optimization scheme that maximizes the smallest
eigenvalue of the Fisher information matrix was applied to the parameter estimation. The
sensitivity coefficients used to construct the information matrix were derived from
temperature histories at sensor locations that, in principle, can be relocated during the
experiment, say through IR thermography, to optimize the estimated parameters. While
IR camera images have been used to determine thermophysical properties in 1-D and 2-D
experiments, the results found thus far have only demonstrated that the technique is
feasible in 3-D, rather than to optimize the experimental design. Phillipi er al. (1995)
used full-field IR images to compute the in-plane diffusivity of an unspecified sample
having a diffusivity of approximately 10'6 mz/s. While a sensitivity analysis showed that

noise in the temperature measurements would not affect identification of the

13

thermophysical parameters, no attempt was made to optimize the experiment to improve

confidence in the results.

A complete analysis of the extended flash method should include optimization of
physical parameters (e.g., location of temperature sensors), as well as experimental
verification of the design of experiment analysis. The research reported here does so by
optimizing the measurements using data available from IR thermography and a variable
power diode laser. The focus here is on optimizing 3-D thermal diffusivity
measurements using D-optimality (Beck and Arnold, 1977) and specifically, optimizing
the sensor position for materials having an out-of—plane diffusivity typical of an elastomer
and in-plane components varying from one-half to ten times the out-of-plane value.
Thermo-optical measurements in atmospheric air are chosen because such boundary
conditions are amenable to biological soft tissues and elastomers, which may outgas in a
vacuum. The unknown and variable boundary conditions add two nuisance parameters
that must be estimated. Experimental demonstration of the optimization technique is
made by positioning sensors for a material having an isotropic response typical of a
metal, in this case, AISI 304 stainless steel plate. Thus, the goal of this study was to
develop a methodology that allows for the optimizing of the temperature measurement

positions post-experiment.

l4

2.1 Methods

The temperature response in a specimen of general anisotropy may be described by the
diffusion equation written as

[ac-35:— = V(kVT) (2.1)

where p is the mass density, c the specific heat, T temperature, t time, and k the thermal
conductivity tensor. The thermal conductivity tensor has nine components, in general,
and is symmetric (i.e., k3 = k32, etc.) such that only six components are independent
(Powers, 2004). The components of k may be considered as being constant over time for
small temperature excursions such as those modeled here. When anisotropy is
considered, the principal axes of k are often assumed to coincide with the axes used to
define the boundary conditions so that only the diagonal terms are nonzero (e.g.,
Mulholland and Gupta, 1977; Doss and Wright, 2000; Venerus et al., 2001). If the
material response is isotropic, the diagonal elements are identical. The thermal
diffusivity a, where a = k/pc, is often more convenient to measure than k because a is

measured during transient, rather than steady state, temperature fields.

Estimating thermophysical properties from temperature measurements requires repeated
solution of the diffusion equation while sequentially reducing the difference between the
modeled and measured results. Here, an alternating direction implicit (ADI) finite
difference solution of equation (2.1) with appropriate boundary conditions is compared to
the bottom face temperature history of the specimen as measured at three points. Figure

2.1 illustrates the geometry of a specimen. Specimens are approximately 40 mm by 40

15

20 mm

 

4\
\

 

 

 

 

 

 

 

 

 

5.3 mm

 

o

 

\/

 

 

 

 

 

 

 

 

A;

Figure 2.1 Diagram of the top surface of the slab specimen illustrating the illumination
area and the specimen geometry. The quarter plane modeled in the parameter estimation
scheme is bounded by the positive x- and y-axes, which themselves are axes of symmetry.
Typically, specimens are between 1.5 mm and 3.5 mm thick. Three sensors are located
on the bottom surface of the specimen along the coordinate axes a distance d apart; the
projections of these sensor positions are shown on the modeled quadrant of the specimen.

mm wide and 3 mm thick. The illuminated area covers a 10.6 mm x 10.6 mm square on
the top face. Due to symmetry along the x- and y-axes, only the area located inside the
positive x- and y-axes is modeled for data reduction. Temperature sensors are located at
three positions: one sensor is located at the intersection of the axes at the center of the
illuminated area, while the other two locations lay a distance d from the center sensor
along the coordinate axes defined by the plane of the specimen. Here, d is to be
optimized for a given specimen geometry and composition, and, in principle, can be

different for the two in-plane directions. For this simulation, however, the x- and y-

16

diffusivities are assumed equal and this distance is assumed to be the same for both in-

plane directions.

The analytical problem models an orthotropic planar specimen that has equal in-plane
dimensions of length 21 and a thickness 8 that is much smaller than the in-plane half-
width 1 (see Figure l). The specimen is heated with a flux qa for a finite pulse time over a
square region on the top (2 = 0) surface. Convection acts over the top and bottom
surfaces, while the side surfaces see small temperature excursions that may be modeled
as adiabatic surfaces. Due to symmetry, only a quarter of the specimen volume shown in
Figure 1 need be modeled. The mathematical formulation of this problem may be stated

as

in - aZT azr azr
pC At—k” Az+k22 A’z-t'kx; /azz (2.2)

mOSxSLOSySLOSzS&t>0

Initial Condition: T(x, y, z, 0) = T... (2.3a)

Boundary Conditions:
—k”3§&giAT-Lg=qomz=dnxosxsnosyshbstsn cam
h,(T-—T,,°)= k33 BT82 atz=0, forlpr SI, l,,_<_y SLOStStp (2.3c)
h, (T - To, ) = k33 3T dz at z = 0, for t > t,, (23d)
3%x=0atx=0andx=l,fort>0 (2-36)
a%y=0aty=0andy=l,fort>0 . (23f)

17

hb(T—T°°)=k33aTaz atz=5,fort>0 (2.3g)

where k,-,- refers to the axial thermal conductivity in the i-direction, T... the ambient
temperature, h, the convection heat transfer coefficient at the top surface, hb the
convection heat transfer coefficient at the bottom surface, 1,, the half-width of the heated
region, and t,, the heating time. During parameter estimation, the thermal diffusivity
components a, (a, = k,, lpc) are estimated instead of the conductivity components kg.
The value of hb, the convection heat transfer coefficient at the bottom surface, is assumed
known, and is based on results using a 1-D conduction analysis. Temperature histories
generated using hb compared favorably to histories generated using a two-layer model
that accounted for heat conduction through the air underneath the specimen, however, the
computational time was reduced significantly by using hb. The ADI model is
unconditionally stable and thus arbitrarily large time steps may be taken during model
solution; however, the accuracy of the stable solution must still be verified. Variation of
the time steps used here resulted in no change in the computed temperature histories
greater than 0.001 K. Studies of grid refinement (in which the grid sizing was changed
from 11 x 11 x 11 to 101x 101x 201) showed that a41 x 41x 41 node mesh accurately
resolved the temperature history of the specimen to within 0.001 K, and the error in the
temperature calculations was within 0.01% based on sz extrapolation. The ADI model
operates within a parameter estimation program that was based on a Levenberg-
Marquardt iterative nonlinear least squares fitting routine (Press et al., 1992). The
parameter estimation program also allowed for spot temperature input from an IR camera

and calculated the confidence region of the estimated thermophysical parameters. The

18

Marquardt program uses a x2 merit function to determine the goodness of fit between the

data and model. The normalized x2 error is defined as

N 2
1
X 2(1); ) = fiZiTj - y 10% i] (2.4)
i=1

where N is the number of data points, p,- the parameters being estimated, T,- the data

points, and yJ-(pi) the calculated result at the point corresponding to the f“ data point.

One criterion commonly used to optimize the design of an experiment, called D-
optimality, minimizes the hypervolume of the confidence interval in the estimated
parameters as expressed by the covariance matrix (Beck and Arnold, 1977). This design
method is well suited for the current analysis, for which the minimization of the
covariance matrix is accomplished by maximizing A, defined as the determinant of the
Fisher information matrix (Beck and Arnold, 1977; Mejias et al., 1999; Emery, 2001).
The Fisher information matrix is composed of normalized sensitivity coefficients Fij (Fij-
= pj(dT,-/Bpj), where T, is the temperature measured at position i and p,— is an estimated
parameter), which indicate the relative change in the temperature field at a point with a
variation in a given estimated parameter. Comparison of different experimental designs

requires defining a scaled Fisher information matrix, which may be written as

1. .. 1 .fim. any. ”ill . f... (.5)
m,n M'tf s-1 {:0 mapm "apt: mer

where M is the number of sensors used (three in this case), If the total time over which

 

 

 

 

data are collected, m and n the indices of the parameters being estimated, T, the

19

temperature at a given sensor and time as calculated by the model using the current
estimate for p,-, and Tm, the highest temperature reached during the experiment. It is
assumed here that the number of temperature measurements taken in the experiments to
be compared is fixed, and that the measurements are sampled at equally spaced intervals
during each experiment. Also assumed are the eight statistical assumptions of Beck and
Arnold (1977), which may be summarized as: errors are additive and uncorrelated with
zero mean and a constant variance, the independent variables are errorless, and there is no
prior information available regarding the estimated parameters. Note that the first two
parenthetical terms inside the integral are Fsm and F5“, which are calculated as part of the
Levenberg-Marquardt estimation scheme, and thus add no penalty to the optimization

calculation.

Although all of the data obtained in an experiment may be used for the optimization,
increasing the information density of the data with respect to the parameters of interest is
highly desirable. That is, there may be times with little change in T at some sensors that
produce little additional information for determining the diffusivity, for example. Huang
(2004) used the ratio of A calculated using all of the estimated parameters divided by A
calculated for the nuisance parameters only, A”, as an optimization criterion (Beck and
Arnold, 1977). As this ratio increases, the uncertainty in the desired parameters (i.e.,
thermal diffusivity) decreases. Because the primary concern is reducing uncertainty in
the desired parameters, this metric provides a clear basis for comparison of experiments
while focusing on the parameters of most importance. Here, optimal sensor positions

will be calculated using both A and N to determine whether optimizing the experiment

20

specifically for the parameters of interest has a noticeable affect on the uncertainty in

these parameters.

To test this procedure for optimizing measurements of the thermal diffusivity tensor,
temperature histories were simulated for an extended flash diffusivity experiment for a
40.0 mm x 40.0 mm x 3.35 mm specimen with an out-of—plane diffusivity approximating
that of PVC (1.08><10'7 m2/s) and transversely isotopic in-plane diffusivities varying from
0.54x10‘7 mZ/s to 10.8x10'7 m2/s, yielding anisotropy ratios 7 of in-plane to out-of-plane
diffusivity of 0.5 < y< 10. Bottom surface temperature histories were calculated at three
points: the center of the specimen, and at two points located on orthogonal axes at a
distance d from the center point, as illustrated in Figure 2.1. The duration of the
simulated experiment and sampling frequency were varied based on the temperature rise
of the in-plane temperature sensors. The inter-sensor distance d was varied from 7.0 mm,
or slightly outside the projection of the illuminated area, to 20.0 mm, at the edge of the
specimen. The inter-sensor position was varied in 1.0 mm increments, except for cases
where the maximums of A and A+ lie between neighboring points, in which case the
sensor position was moved to the 0.5 mm increment bisecting these two locations to
improve interpolation of the optimal sensor location. The heating duration and power
were set such that the maximum temperature of the specimen was 100 °C. While this
temperature is higher than the glass transition temperature of PVC (T8 = 81 °C, CRC
Handbook of Chemistry and Physics, 1991), it is less than the 105 °C regularly
experienced by some rigid PVC tubing regularly used to carry electrical cables

(Matthews, 1996). Moreover, the temperature excursion in the model exceeds the 80 °C

21

level for no more than 0.4 s for a 1.0 s heating time and 5.8 s for a 10.0 s heating time,
and the volume over which T> 80 °C comprises less than 3% of the total specimen
volume. Random noise was added to each computed temperature history to simulate
measurements such as those acquired with an IR camera. The noise had a Gaussian
distribution with zero mean and a standard deviation of o = 0.01 K, which is on the order
of the noise resulting from IR thermography. Additional calculations were performed
using temperature histories that had added noise with standard deviations of 0.1 K and
1.0 K to determine the effect of the noise level on the data reductions. Once the
simulated temperature measurement histories were computed, each served as data for the
parameter estimation program, and the three desired parameters (ax, cry, and a2) and two
nuisance parameters (the incident heat flux and convection coefficient on the top face)
were estimated. The output of the reduction program consisted of the five estimated
parameters, the normalized x2 error for the model fit, A and N, and scaled sensitivity

coefficients and residual plots for each of the three sensor locations.

In order to validate the optimization routine, the optimization was repeated for the sensor
positioning on a simulated specimen having thermal properties similar to those of steel,
and these results were compared to extended flash diffusivity tests conducted on an A181
304 stainless steel plate. This material was chosen for verification because of its
isotropic nature and well-characterized thermal diffusivity (CL = 3.95 x 10'6 m2/s;
Incropera and Dewitt, 1990). The 70 mm x 70 mm x 1.50 mm plate was sprayed with a
thin layer of flat black stove paint to increase absorption of the laser heat pulse on the top

surface of the specimen as well as to decrease the reflectivity and enhance the IR

22

thermographic response on the bottom surface. Four extended flash diffusivity tests were
conducted on the plate using heating parameters developed for steel in prior studies
(LeGall, 2001). Heating was performed using a 16W diode laser with a beam that was
expanded, collimated and then truncated with an aperture plate to achieve the desired
heating area. Four tests were performed, each with a heating time of t,, = 1.0 s. Full-field
data were sampled at the maximum frame rate of 60 Hz for 30 5. Temperature histories
were retrieved from the full-field data at three points on the bottom of the specimen
surface: under the center of the illuminated region, and two spots located (I from that
center spot, as shown in Figure 2.1. Two values of d were chosen to illustrate the relative
ability to estimate parameters at locations both close to and far from the optimal point as
determined by the prior optimization analysis. Temperatures for each sensor location
were averaged over a 3-pixel by 3-pixel area that was centered over the pixel of interest,
in order to reduce random noise. Five parameters were estimated for each test: the three
components of the thermal diffusivity tensor (ax, 0ty and a2), a scaled heat flux
coefficient, q/pc, and a scaled top surface convection coefficient h/pc. The average
values and standard deviations of the out-of—plane diffusivity, a2, and the in-plane
diffusivities, 0tx and (1,, were estimated for each sensor position and compared to the
values of A and AI calculated from the simulated runs and the measurements to ascertain

whether a correlation exists between the accuracy and variation in the values of aand the

magnitudes of A and A" for each set of sensor positions.

23

2.2 Results

2.2.1 Measurement Duration Figure 2.2 shows simulated temperature histories and
corresponding inverse solutions for four cases of an isotropic material with or = 1.1 ><10'7
m2/s, approximating PVC. Panel (a) shows temperature histories for t,, = 1.0 s and two x-
axis sensor positions, d = 7.0 mm and d = 20.0 mm. Panel (b) shows similar temperature
histories for to = 10 s and sensor positions of d = 7.0 mm and d = 20.0 mm. One second
corresponds to a Fourier number based on the specimen thickness of Foz = (020/52 = 9.62
x 10'3 and that based on the specimen length of Fox = (0t,,t)/l2 = 6.75 x 10“, where 8 is the
specimen thickness and l the length of the specimen side. Because the materials modeled
in this study are in general transversely isotropic, the y-axis temperature histories are the
same as the representative x-axis histories, within the magnitude of the added noise, and
are omitted here. Calculations for t,, = 0.1 s and t,, = 10 s were conducted for the isotropic
specimen to give insight into the relative magnitude of information derived from each set
of tests. The inverse solution fits the simulated experimental results well, as shown by
the residuals, defined as the difference between the temperature histories of the simulated
experiment and the inverse solution (Figure 2.3). Ideally, the residuals should oscillate
about zero with a deviate on the order of the noise exhibited in the temperature

measurements.

Comparison of A and A+ for different experiments requires each test to be designed so
that the same number of data points are collected for each trial (Mejias et al., 1999).

Such design allows comparison of experimental design without skewing the results by

24

 

9 (K)

 

 

 

 

 

 

-1 . , . . .
0 50 100 150 200
time (s)

 

 

 

 

 

 

-1 . , . . .
50 100 150 200

time (s)

Figure 2.2 Simulated temperature histories and inverse model reductions for an isotropic
Specimen at sensor locations of d = 0, 7, and 20 mm: (a) 1‘p = 1 s, (b) t,, = 10 s. O is
defined as T-T... where To. is the ambient temperature.

25

 

0.04

0.02

Residual (K)
O
O
O

 

-0.02

 

 

-o.o4 . . t . t t
0 50 100 150 200
time (s)
0.04
(b)
0.02

Residual (K)
O
8

 

-0.02

 

-0.04 1 . I i i
O 50 100 150 200
time (s)

Figure 2.3 Representative residuals, defined as the difference between the modeled and
measured temperature histories, for an isotropic specimen: (a) center temperature sensor,
(13) off-axis temperature sensor.

26

using different quantities of data. A common sampling technique to achieve equal
numbers of data points is to use the same acquisition rate for the same duration over the
experiment. Changing material or experimental protocol, however, may result in these
experiments collecting large amounts of relatively unimportant data. The temperature of
a thick, low diffusivity specimen, for example, may require hundreds of seconds to reach
ambient temperature following a step heat input while a metal specimen may reach
ambient temperature after only tens of seconds. Similarly, adjustments in positioning the
sensors may yield temperature histories with different times of maximum signal.
Consider the temperature response in Figure 2.2a. At the 7.0 mm sensor, the temperature
begins to decrease after 65 s, but the temperature at the 20.0 mm sensor has barely risen
after 200 5, suggesting that a faster sampling rate and shorter duration are required for the
7 mm sensor, as compared with the 20 mm sensor, to retrieve similar amounts of useful
information. Instead of fixing the sampling rate and duration, it may be advantageous to
fix the number of data points sampled during the temperature transient and vary the
sampling rate and duration based on some other criterion. A series of computational
experiments was run for an isotropic specimen to investigate the quality of information
taken from such a test. Each test contained 2000 data points at each of three sensors over
the duration of the experiment (for a total of 6000 points). The sampling rate and
duration were varied to yield 6000 points from a minimum duration of about 20 s to a
maximum of 250 5, thus capturing the temperature histories from the initial temperature
rise on the bottom center of the specimen to the cooling of the bottom specimen surface.
Three inter-sensor distances were chosen to give insight into the effects of changing the

sensor location on the optimum sampling duration. Figure 2.4a shows the in-plane

27

temperature response of an isotropic specimen for d = 12 mm and the corresponding
sensitivity coefficients. Figure 2.4b shows the values of A and A+ as a function of elapsed
time for the same sensor location. Note that the maximum values of A and A+ occur at
approximately the same time as the maximum value of the scaled heat flux sensitivity
coefficient qu, as shown in Figure 2.5. For tests conducted using small d, these maxima
occur slightly after the in-plane sensors reach their maximum temperature, but as d
increases, the maximum of the in-plane temperature response occurs at approximately the
same time at which A, N, and qu reach their maximum values. These results suggest
that sampling the temperature histories until each of the in-plane temperature sensors has
reached its maximum value should yield the most robust results for each test. Extended
flash diffusivity measurements on PVC and steel by Doss (1998) and LeGall (2001)
corroborated these results, as the variability in the estimated diffusivities was minimized
when the measurements were conducted until maximum rise of the in-plane temperature

sensors.

2.2.2 Anisotropic Diffusivity This approach has been extended to materials with an
orthotropic thermal response. Simulations without added noise were used to identify the
time of the maximum in-plane temperature excess for 0.5 S y S 10 and 6 mm S d S 20
mm. Once these times were identified, data were simulated again by adding noise to
temperatures calculated using the parameters noted above. The sampling rate and
duration were adjusted to yield 6000 data points (2000 points from each of the three
sensor locations) over the calculated time of maximum temperature. These data were

then used as input for the parameter estimation scheme, which again estimated five

28

 

 

 

 

 

 

 

 

 

1.0+ (b) 2‘ 0

OD

OD
O

 

 

0.0- u‘

 

0 50 100 150 200 250
time (s)

Figure 2.4 Panel (a) shows the in-plane temperature response and Hg for an isotropic
specimen where d = 12 mm. The dots represent the in-plane temperature data, the solid
line the CL, sensitivity coefficient, the dotted line the heat flux sensitivity coefficient, the
short dashed line the (1.1 sensitivity coefficient, the long dashed line the or, sensitivity
coefficient, and the dot-dash line the convection sensitivity coefficient. Panel (b) shows A

(circles) and A” (triangles) as a function of elapsed time for the same specimen and
sensor location.

29

 

-7 ' . 4
' ' 0‘ t ' 00
.v
“ J
.a - -.. 0.
fl 1.! ‘

¢ .

’v‘.
D.

l

A

V

J
"\.- ' '. .
v ' . I'I -
. O o
' - . 1 '
.
-' . -.I . .‘ - p. n, t
"' .. .' .. -~.--.~ ~.« I . y—
n - p . 5 ‘. - a, - ,‘
v~\. I fit . - o
_-.- 0 f . '
w
... Q ‘ L. .
o ! '

\ 4 .
't ‘ ' I
‘ .,_, .J
c I <_“
n" .' P
'.‘«:‘<
.“
'3"
"3:
9 — ‘Kfi. " p— o 2
- - __
. . o
.

log A+
9 (K)

-1 0 ‘ 39:" “ 0.1

-11 4E7 ' 0-0

 

 

 

-12 ‘ I ' l ' l ’ ' T '01
0 50 100 150 200 250
time (s)

Figure 2.5 In-plane temperature response (circles) and log A+ (triangles) as functions of
time for an isotropic specimen with d = 12 mm.

30

parameters and calculated A and A“. Figure 2.6 illustrates the typical variation of A and
A+ for changing d, with the maximum for each A and A+ lying at d = 8 mm. For each y,
the out-of—plane diffusivity and the average of the two in-plane diffusivities were
resolved to within less than 1%. Figure 2.7 shows how the optimal d varies with y as
described by the maximum A and N. These maxima occur for sensor positions near the
illuminated region for small 7, and closer to the edge of the specimen as 7 increases.
Note also that the maximum values of A and A” occur at the same location for y < 6, but
separate for larger values of 7. For 7 = 9 and 10, the maximum values of A and A+ occur
at d = 20 mm, the edge of the specimen, suggesting that the specimen width must be
increased to accommodate materials with y > 8. Due to the ambiguity in the meaning of

these results, only results for y S 8 will be discussed further.

2.2.3 Experimental Validation Extended flash diffusivity tests were conducted on
A181 304 stainless steel to demonstrate the accuracy gained in the estimated parameters
when using optimized sensor positioning. The analytical optimization routine detailed
previously was repeated for a simulated 1.50 mm specimen having an isotropic thermal
diffusivity (or = 3.95 x 10'6 mzls) and heat flux and convection boundary conditions
similar to those in prior flash diffusivity tests (LeGall, 2001). As shown in Figure 2.8,
the simulated experiments indicated that the optimal sensor position as quantified by both
A and A+ should be as far from the illuminated region as possible, suggesting that d
should be set to 20 mm for a 40 mm x 40 mm specimen. In principle, such an
experiment might yield the results with highest confidence, except that the assumption of

the specimen edge being adiabatic is violated, as there is a small temperature gradient at

31

 

 

 

 

0 o
. o
-4* . ‘ --8
‘ o

o

<’ “<1

8’51 ‘ ”-93
A o

'5‘ . --10

A o
o
'7 I'l'I'TTFTI'I'I‘1*I'I'I'I*T‘-11
6 7 8 9101112131415161718192021

d(mm)

Figure 2.6 Variation of A and A+ with respect to d for an isotropic specimen with 1,, = 1.0
8. Calculated values of A are represented by circles, while calculated values of A“ are
represented by triangles.

32

 

21
19‘
17 ‘

15‘

 

 

 

Figure 2.7 Location of the maximum A (circles) and A" (triangles) as a function of y.

the specimen edge. Additionally, prior measurements with the IR camera suggest that the
temperature rise at the off-axes sensors is on the order of the experimental noise for d >
15 mm and a heating time 1,, = 1.3 s. In light of this, sensor positioning was limited to the
area bounded on the interior by the shadow of the heated area and the exterior by roughly
the outer third of the specimen half-width and, thus, 6 mm < d < 14 mm; specifically, d =
8.45 mm and d = 13.75 mm were chosen for the IR measurements. These values of d
correspond to distances of 40 and 65 pixels, respectively, from the central illuminated
pixel in the IR camera field of view. Results of these measurements are shown in Table

2.1, Table 2.2 and Table 2.3. Note that the average of the estimated 0t" and ory using d =

33

8.45 mm differs from the reference value by almost 5%, but the difference between the
reference and average estimated in-plane value using (I = 13.75 mm drops to 0.5%. The
error in the estimates of 0tz are higher in magnitude, over 11% for d = 8.45 mm, but about
5% for d = 13.75 mm. Table 2.3 lists the experimentally obtained values of A and AJ, for
the steel plate tests, which are, in general, lower than the values calculated from the
analytical simulations. The estimated values of the nuisance parameters (heat flux and
convection coefficient) were slightly lower than those used for the simulation, resulting
in lower calculated values of Fij, A and N. For tests conducted at the same d, the average
values of A and A+ are different for different heating durations. A longer heating duration
creates higher temperature excursions in the specimen and, therefore, larger scaled
sensitivity coefficients, increasing the value of A. Note that the average values of both A
and A+ are about three to four times larger for the d = 13.75 mm tests, where the
estimated diffusivity values are closer to the reference values, as compared to those

values obtained for d = 8.45 mm.

34

 

 

 

 

'4 -8.5
,
'5_ ‘ A
. A . O
. ' --9.3
-6“ ‘ .
A ‘ . .
<1 '7 . +<,
U) --
2 A . 10.0 .8)
-8“
C
.9n
”-10.8
40* '
'11'1'I'I'r'l'i'r*r't'1't'tvrtrv-11.5
6 7 8 9 1O 11 12 13 14 15 16 17 18 19 20 21
d(mm)

Figure 2.8 Variation of A and A+ with respect to d for an isotropic steel specimen with t,,
= 1.0 5. Calculated values of A are represented by circles, while calculated values of A+
are represented by triangles.

35

Table 2.1. Estimated directional thermal diffusivity values and associated standard
deviations for a 1.50 mm stainless steel plate (average of four tests).

(I

 

 

 

 

 

 

 

avera e orz on average (itx ow, average 0ty on,
(mm) (mm ls) (mmzls) (mm2/s) (mm2/s) (mmZ/s) (mm ls)
8.45 l 3.50 ‘ 0.13 l 3.74 i 0.07 , 3.78 l 0.05
13.75 3.75 0.14 4.01 0.06 3.93 0.04

Table 2.2. Average deviation of the out-of-plane and in-plane diffusivities from the
reference value of or = 3.95 x 10’6 mzls for a 1.50 mm stainless steel plate.

 

 

 

61 average or, deviation average 01,, 0ty deviation from
(mm) from reference (%) reference (%)
8.45 I -11.4 ’ —4.8
13.75 -5.1 0.5

Table 2.3. Values of A and A+ computed from measurements of the 1.50 mm steel plate

d (mm) I

 

average A I average A+
8.45 3.02 x 10‘" 2.43 x 10‘l2
13.75 1.30 x 10'9 8.03 x 10'12

36

2.3 Discussion

The estimated values of thermal diffusivities in the analytical portion of the study all fell
within 0.2% of the input diffusivities for 1 S 7S 8 with added noise having a standard
deviation of 0.01 K. The uncertainty in the estimated diffusivities reported in the
literature for multiaxial measurements is approximately an order of magnitude larger than
the uncertainty calculated here. For example, measurements made by Doss and Wright
(2000) on A181 304 stainless steel exhibited a standard deviation of 3.5% of the reference
diffusivity value in the out-of-plane direction and an uncertainty of 5.6% over the in-
plane directions, while 2-D measurements on PVC by Maillet er al. (1990) showed a
variability in the estimated diffusivity value of 3% between tests in the axial (out-of-
plane) direction and 5% in the radial (in-plane) direction. Doss and Wright reported
noise on the order of 0.1 K, which may result in the greater uncertainty. The higher
accuracy of the computational experiments may be expected, because the same analytical
model that generated the temperature histories, after addition of noise, was also used in
the parameter estimation, and thus the simulated experimental boundary conditions are
perfectly modeled. Moreover, the noise in the simulations conforms to the eight standard
assumptions for parameter estimation (Beck and Arnold, 1977). For 7 = 0.5, accuracy
decreased dramatically for d > 13 mm, with the error in the estimated values of the
thermal diffusivities exceeding 4% or more for the out-of—plane direction and 70% for the
in-plane directions. For these cases, the in-plane temperature rise was less than 0.07 K,

leading to low signal-to-noise ratios and sensitivity coefficients of low magnitude. The

37

lack of a strong signal results in poor quality information at those sensor locations and

thus, the inability to estimate the parameters properly from such data.

The information obtained from each experimental simulation, as quantified by A for the
optimal sensor position and heating time, varied from a maximum of 1.35x10'3 for y =
0.5 to 1.00><10'6 for y = 8. For A”, a maximum value of 1.58x10'8 was calculated for y =

1 and a minimum value of 2.83><10"0 for y: 8. Table 2.4 shows the maximum values of
A and N along with the corresponding (1 as a function of 7. Note that both optimality
criterion yield identical d for 1 < y < 2, a range typically seen in deformed or machined
elastomers (e.g., Doss and Wright, 2000; Venerus et al., 2001). One might reason that
the closer the sensors are to the heated area, the higher the temperature signal and,
presumably, the more information is obtained from the experiment. Results based on A
and A’( indicate that this is not necessarily the case, however, as the optimal sensor
position may be located far from the illuminated area, as shown in Figure 2.7 for
elastomeric specimens having a large 7. For 7 < 6, the optimum sensor position lies at
about the same location, either 8 or 9 mm, for either criterion. For 6 < y < 8, however,
the optimal position moves further from the center, and this location differs depending on
whether one bases optimization on the maximum A or A+ (see Table 2.4). While prior
studies have looked at the optimization of sensor positioning (e.g., Maillet et al., 1990;
Mejias et al., 1999; Mzali et al., 2004), the optimal position was not examined as a
function of anisotropy ratio as it is here. Optimizing the measurements as a function of
varying anisotropy is important because the thermophysical properties of elastomers have

been shown to vary with mechanical loading (e.g. Broerman et al., 1999). Since

38

Table 2.4 Optimum sensor location based on A and A“ for various y in a phantom
elastomer with an out-of—plane diffusivity of l.08><10'7 m2/s and 1,, = 1.0 s.

 

 

y optimum A d for optimum A optimum A+ d for optimum A+
(mm) (mm)
0.5 1.35 x 10'3 8 1.32 x 10'8 7
1.0 7.56 x 10“ 8 1.58 x 10'8 8
2.0 1.49 x 10“ 9 7.63 x 109 9
3.0 5.00 x 10'5 9 4.50 x 109 9
4.0 1.56 x 105 9 2.14 x 10° 9
5.0 6.38 x 10*5 10 1.27 x 109 10
6.0 2.79 x 10'6 12 6.78 x 10'10 10
7.0 1.43 x 1045 13 4.19 x 10'10 11
8.0 1.00 x 106 16 2.83 x 10'l0 12

 

 

 

material thermophysical properties may be unknown prior to estimation, especially if the
material is deformable and affected by loading during the experiment, repositioning of
the experimental temperature sensors based on the measured properties is necessary to
yield the most accurate results. Methods of temperature measurement that allow
flexibility in sensor positioning post-data acquisition enable an iterative data analysis that
will yield better estimates of thermal pr0perties than those afforded by fixed-position

methods of temperature measurement, such as thermocouple probes.

Variation of the simulated heating time seemed to have little effect on the location of the
optimum sensor position. Simulations using 1,, = 0.1 s were inconclusive, as the
simulations for d > 1 1.0 mm did not converge for that heating time due to the low
temperature signal resulting from the imposed 100 °C temperature limit on the top face of
the simulated specimen. Table 2.5 shows the relationship between t,, and the maximum
values of A and A”, as well as the maximum error in the estimated diffusivity values for

each heating time. While the information available to the reduction program for 1,, = 0.1 s

39

was enough to resolve the diffusivities to within 1% for d S 11 mm, the differences in
magnitude of A and A” near the optimum were very small (less than 3%), leaving the
optimal d unclear, although it appeared to be between 8 mm and 9 mm. Results for 1,, =
10.0 5 indicated an optimal sensor position of 8.0 mm based on A”, similar to the 1.0 s
heating time tests, but the optimal sensor position using A shifted to 9.0 mm. As the
heating time increased, the optimum values of A and A” increased by orders of
magnitude, suggesting that the accuracy of the estimated parameters would improve with
increasing heating time. Prior measurements have shown that extended heating (e.g., FoZ
> 0.2 corresponding here to t> 20 5), even at low power, can complicate the modeling of
the convection on the specimen surface, as well as cause large temperature gradients near
the specimen edges, both of which depart from the modeled behavior (data not shown).

As a result, simulations having longer heating times were not investigated, despite A and

Table 2.5 Relation between t,,, optimum A and N, and variation in or.

 

’1) (S) optimum A optimum A+ maximum error in 0t (%)
0.1 3.31 x 10‘8 2.73 x 10'” 0.5
1.0 7.56 x 10“ 1.58 x 10'8 0.1

10.0 6.50 x 100 1.43 x 10‘5 < 0.1

 

 

 

A+ both increasing as 1,, increased. Note that the maximum variation in the estimated
diffusivity values in Table 2.5 lies below one percent, even for the cases with short
heating times. Given the accuracy possible with the 1,, = 0.1 5 case, one might conclude
that the shortest heating time that gives a clear signal may be acceptable for running such
transient thermophysical property measurements. Experiments using essentially

instantaneous flash times (Degiovanni, 1977) may require such a large heating power (to

40

obtain a good temperature signal at all sensor positions) that the top face of some
specimens may ablate or undergo a phase change, which would be undesirable. Thus, if
the total energy input to a specimen is to be held constant, the longest possible heating
time should be used that minimizes edge effects or complicated convection boundary

conditions.

The level of noise present in the isotropic data affected the accuracy of the estimates for
the diffusivities, but did not affect the optimum value of d. Near the optimum d (:2 mm
from the optimum location), the accuracy of the estimated diffusivities was
approximately the same for both the out-of-plane and in-plane directions, but the error in
the estimated out-of—plane diffusivity was approximately three times lower than the
associated error in the in-plane directions when the sensor position was far from the
optimum d (see Table 2.6). Note that, for the case of maximum noise, when 0' was equal
to 1.0 K, the variation in the out-of-plane and in-plane diffusivities was at most i10.6%
when d was close to the optimum, but this variation increased to i27.7% for the in-plane
diffusivity far from the optimum configuration. While the magnitude of A at the
optimum location decreased by 47% as the standard deviation of the noise increased from
0.01K to 1.0K, the magnitude of N at the same points remained nearly constant. Because
AI relates the ratio of the overall information divided by the nuisance information, a
decrease in the overall information A, with A” remaining constant indicates that an
increase in noise level has a negative effect primarily on the quality of the information in
the nuisance parameters. While attempts should always be made to minimize noise and

improve the signal-to-noise ratio when conducting flash diffusivity experiments, these

41

Table 2.6 Relationship between noise level in measurements and the corresponding
accuracy in the thermal properties of interest for an isotropic specimen with t,, = 1.0 s.

 

 

 

 

 

0' (K) maximum variation for (itz maximum variation for 01., far optimum value
near optimum d (%) from optimum d (%) of A
0.01 < 0.1 0.1 7.56 x 10'4
0.1 0.5 0.5 7.00 x 10“
1.0 4.8 6.1 3.97 x 104
c (K) maximum variation for maximum variation for Other, optimum value
abut, near optimum d (%) far from optimum d (%) of A+
0.01 < 0.1 0.5 1.58 x 10'8
0.1 0.6 3.5 1.59 x 108
1.0 10.6 27.7 1.58 x 10’8

 

 

 

results indicate that extreme efforts to reduce noise in the experiments may not result in
an improvement in the accuracy of the estimates of the thermal diffusivities, since the
noise primarily affects the estimation of the nuisance parameters. However, the use of IR
thermography to measure temperature has reduced the level of noise by an order of
magnitude as compared to E-type thermocouples for measurements in our laboratory, and
this noise reduction corresponded to a reduction in the number of iterations needed to
converge on the estimated parameters. Thus, the benefits of noise reduction may extend
beyond the accuracy of the parameters. While the analytical results for y < 6 suggested
that either A or A” could be equally useful in determining the optimal d, the reduction in
computing time brought on by the improved signal-to-noise ratio suggest that the relative
magnitude of experimental noise affects the certainty in the parameters of interest via a

decrease of information in the nuisance parameters, and as such, A may be more useful

than A+ for experimental optimization.

42

Experimental validation of the optimization procedure was performed using a stainless
steel plate. Analytical results showed that more information was obtained about the
parameters of interest (as quantified by A and A‘“) at the far temperature sensor location
despite the overall lower temperature signal at that position. The increase in information
obtained at the optimal sensor position resulted in estimated diffusivity values that were
closer to the reference diffusivity of A151 304 stainless steel than diffusivities estimated
using the lower information available at the non-optimal sensor location. The standard
deviation of the estimated in-plane diffusivities compared favorably to earlier
experimental results by Doss and Wright (1.3% of the reference value for these
measurements versus 5.6% for Doss and Wright), while the standard deviation in the out-
of-plane diffusivity 0tz was unchanged (3.5% in both sets of experiments). However, the
estimated value of (1,, reported here was 5.1% lower than the reference value, while the
value reported by Doss and Wright was only 1.8% lower than the reference value. This
discrepancy in 0t; is likely due to the 60 Hz frame rate of the IR camera. The camera
frame rate is constant for a given test and does not synchronize to accuracy greater than
60 Hz with the heat pulse trigger, and as a result there is some scatter in the estimated
out-of-plane diffusivities over successive tests. The temperature response at the off-axis
sensor locations is slower, and as such, a temporal error amounting to a fraction of a time
step does not affect the estimated in-plane diffusivity values. This temporal effect is
magnified for materials of high thermal diffusivity, such as steel. When the steel data
temperature histories were adjusted by adding or subtracting a single time step of 1/60 s,

the estimated value of orz varied by 17.5%. When similar adjustments were made to the

temperature history of a material having a smaller thermal diffusivity (PVC sheet; at z 1.3

43

x 10’7 mZ/s), the estimated value of 0tz varied by less than 1% despite the use of a slower
frame rate of 15 Hz to capture the PVC data. Because this timing error is random, it is
expected that the average value of (1tz will approach the reference value as the number of
measurements in the sample is increased. Still, the results gained from the optimal sensor
location reduced uncertainty in the experimentally derived values of diffusivity, and as
such, shows that this method holds promise for improving future extended flash

diffusivity measurements.

2.4 Closure

Estimation of the axial components of the diffusivity tensor was possible to within a
0.2% error for all experiments featuring a 1.0 s heating time, an CL, = 1.08x10'7 m2/s, and
l < "y < 10. The optimum d was determined via two D-optimality parameters: A, which
gives a measure of the total information available from an experiment, and A+ which is
the ratio of A as defined above divided by A for just the nuisance parameters (here, heat
flux and convection). This distance d was shown to be a function of y for the analytical
experiments simulated here. Increasing the pulse time results in both a higher level of
total information and a higher level of information concerning the parameters of interest,
but physical complications to the experiment may outweigh the slight improvement that
results in the accuracy of the estimated parameters. Experimental measurements made on
a stainless steel plate showed that measurements made closer to the optimal sensor
positions as suggested by this analysis lead to estimated parameters that are closer to the
expected reference values. These results support the use of IR thermography, which
allows for more flexible sensor positioning and lower measurement noise than other
fixed-position methods of measurement for transient extended flash thermal diffusivity
measurements. Overall, with this method showing promise in optimizing experiments on
specimens having transverse isotropy, a next step would be to investigate materials

possessing thermal orthotropy.

45

CHAPTER 3

ON MEASURING THE THERMAL DIFFUSIVITY TENSOR

WHEN PRINCIPAL AXES ARE UNKNOWN A PRIOR!

3.0 Motivation

When modeling the anisotropic thermal response of materials, the material may be
rotated such that the physical axes coincide with the principal axes of the thermal
diffusivity tensor, resulting in thermal orthotropy. This may be impractical for a
specimen of arbitrary shape. During measurements of such a tensor, however, the
principal axes may be unknown, requiring a method to determine principal values and the
orientation of the principal directions while simultaneously measuring the diffusivity.
Thermal anisotropy may result from the microstructure of the material (e.g., crystallinity,
grain boundaries, molecular alignment in polymers) or the macrostructure (e.g.,
composite fibers, wood grains). Deformation and manufacturing may induce anisotropy.
Calendaring of PVC, for example, can produce differences of 30% in the components of
thermal diffusivity of an apparently isotropic material (Doss and Wright, 2000).
Measuring the thermal diffusivity tensor may thus provide a nondestructive test method
for the material response to a mechanical state if the thermophysical properties can be
correlated to deformation (LeGall, 2000; Choy et al., 1978) or stress (Venerus et al.,

1999).

46

The value of each component of the thermal diffusivity tensor depends on the reference
orientation from which the specimen is viewed. The temperature response in a specimen
of general anisotropy may be described by the diffusion equation, written as

pc%§=V(kVT) (3.1)

where p is the mass density, c the specific heat, T temperature, t time, and k the second-
order thermal conductivity tensor. The thermal conductivity tensor, in general, has nine
components:

Fkr 1 kn k13 1

k= [<21 1‘22 k23 (3.2)
H31 1‘32 1‘33

 

 

k is a symmetric tensor with only six independent components, as required by Onsager’s
theory of reciprocity, and thus kij = kji (Onsager, 1931; Carslaw and Jaeger, 1959; Powers,
2004). For small temperature excursions, the elements of k may be considered constant
over time. The thermal diffusivity a, where a = k/pc, is often more convenient to

measure than k because a is measured during transient, rather than steady state,

temperature fields. The diagonal terms of 11(01“, an, (133) are associated with axial
conduction in Cartesian coordinates, while the off-diagonal terms, here called the
shearing diffusivities, define heat conduction along the lines that bisect the axes indicated
in the mixed derivative dleaiaj associated with each cg,- term. The effect of the shearing
diffusivity acting in the 1-2 plane is illustrated in Figure 3.1. Figure 3.1a shows
temperature contours on the bottom surface of a specimen having equal axial diffusivities
in the 1- and 2-directions with no accompanying shearing diffusivity, resulting in

transverse isotropy, while Figure 3.1b shows the same contours for a specimen having the
47

same axial diffusivities as before but now also possessing a shearing diffusivity that is
about 18% of the axial diffusivity value. Note how the addition of the shearing
diffusivity causes elliptical-shaped temperature contours with major and minor axes
oriented towards the lines y = x and y = -x, as opposed to the circular contours that would
result from transverse isotropy (as shown in Figure 3.1a). The change in isotherm shape
due to the effects of the shearing diffusivity is analogous to the distortion experienced by
a material element that is deformed under a shearing strain in solid mechanics, hence the
term “shearing” diffusivity. If the principal axes of thermal diffusivity coincide with the
physical axes of the experiment, the off-diagonal shearing terms have zero magnitude,
resulting in material orthotropy (Figure 3.1c). Rotation of the case presented in Figure
3.1b clockwise by 45° results in the same contours shown in Figure 3.1c. The
temperature response depicted in Figure 3.1b may result from an orthotropic material
having principal axes that are rotated an angle B from the physical axes of the

experiment.

Several studies have presented analytical models of orthotropic materials useful for
thermophysical property estimation. Mulholland and Gupta (1977) discuss the 3-D
solution of the diffusion equation for an anisotropic solid with boundary conditions that
are difficult to solve in an orthonormal coordinate frame. The model is rotated such that
the principal and physical axes coincide and then the coordinate lengths are scaled by a
factor that renders the resulting equation isotropic. The transformation of an orthotropic
material response into an isotropic response via appropriate scaling of the coordinate axes

is discussed further by Ozisik (1993). Chang and Tsou (1977a, 1977b) developed an

48

20

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

20
(a) , (b)
10‘ 10*
E ’5‘
e 0* e 0:
> >
-1o< -10‘
-20 . . . _20 . . r
-20 -10 O 10 20 -20 -1O 0 1O 20
X (mm) 20 x(mm)
(C)
10‘
’e‘
e 0‘
>
-10'
-20 . . .
-20 -10 0 10 20
x(mm)

Figure 3.1 Isotherms on the bottom surface of a specimen at t = 30 8 following 1 s
heating. (a) transversely isotropic: or” = 0122 = 0.16 mmzls, 01.2 = 0, (b) anisotropic: 01“ =
0122 = 0.16 mmz/s, (1.12 = 0.0282 mmzls, (c) orthotropic: or“ = 0.188 mmzls, 022 = 0.132

mmzls, 002 = 0. Note that Figure 3.1c is congruent to Figure 3.1b if rotated 45°
counterclockwise.

algebraic transformation using conformal mapping that reduces a thermally anisotropic

problem to an isotropic one, but this technique was restricted to problems formulated in

cylindrical coordinates having homogeneous boundary conditions. Poon et al. (1979)

extended the conformal mapping technique to specific 2-D and 3-D cases in both

49

rectangular and cylindrical coordinates with homogeneous boundary conditions. Huang
and Chang (1984) extended this conformal mapping technique to 2-D geometries having
mixed boundary conditions (i.e., flux at one boundary, fixed temperature at another).
Haji-Sheikh and Beck (2002) developed an analytical solution for the temperature field in
3-D, multi-layer bodies where the principal thermal and physical axes coincide, leading
to thermal orthotropy, and thus, three nonzero terms in the 01 tensor for each material
layer. Aviles-Ramos (2002) and Mejias et al. (1999; 2000; 2003) have investigated the
optimization of sensor position and experimental boundary conditions for orthotropic
specimens. Pron and Bissieux (2004) have developed a 3-D model of heat diffusion in
orthotropic media for investigation of stress-induced orthotropy. Mzali et al. (2004)
analytically investigated the effect of sensor positioning on the estimated diffusivities
involving an orthotropic cylindrical specimen. Recent analytical treatments of property
measurements on anisotropic thin films (Bennett, 2004) and a 2-D cylindrical slab
(Milosevic and Raynaud, 2004) have assumed that there is no shearing conductivity that
acts between the axial and radial directions, and furthermore, that the radial conductivity
is isotropic with respect to in-plane orientation. Powers (2004) lists necessary conditions
on the conductivity tensor that allow only physically meaningful anisotropic problems;
namely, the conductivity (or diffusivity) must be positive semi-definite and that
Onsager’s reciprocity theory must hold (00,- = 09,). Aside from a treatment of a steady-
state 2-D solution for slab-like multilayer samples (Hsieh and Ma, 2002; Ma and Chang,
2004), little work exists in the recent literature that investigates the response of materials
oriented in other than the principal directions, which yields a fully populated diffusivity

tensor.

50

Prior experimental techniques used to measure the thermal diffusivity tensor have
assumed that the material response is orthotropic, resulting in nonzero values for only the
diagonal terms in 0. Early extensions of the flash method that looked at both an axial
and radial component of thermal diffusivity assumed no coupling between these two
terms (Donaldson and Taylor, 1975; Chu et al., 1980) and thus no shearing diffusivity
need be considered. Doss and Wright (2000) extended the 1-D flash method of Parker et
al. (1961) to three directions through estimation of the two principal in-plane components
of or via the temperature history of a specimen at three distinct points. It was assumed
that the material either exhibited no shearing diffusivity or that the principal axes of a
were aligned with the physical axes of the experiment, resulting in thermal orthotropy.
Graham et al. (1999, 1999b) suggested an experimental format to measure four of the six
independent components of the thermal diffusivity tensor in a slab specimen. Their
formulation assumed that the out-of-plane direction is a principal direction, and thus there
were no shearing terms involving the out-of—plane axis. The in-plane orientation was
assumed unknown, and it was suggested that ratios of temperature histories obtained
from the top and bottom surface of an anisotropic specimen may be used to ascertain both
the three principal diffusivities and the angle B that the principal in-plane axes make with
the physical axes of the specimen. No data were presented, however, to illustrate the use
of this experimental design for thermal property estimation. Also, while Graham et al.
did estimate the sensor location in each analytical experiment, the experimental sensor
positioning itself was not optimized. Because the assumption of thermal orthotropy

seems prevalent in recent experimental studies, quantification of the error inherent in

51

assuming orthotropy for parameter estimation studies in a material exhibiting randomly
oriented anisotropy would be of value. Such a study could establish conditions for which

orthotropy could be assumed in a material exhibiting a fully populated 0t tensor.

To investigate the effects of randomly oriented anisotropy on the temperature response of
a specimen, computational models need to be written that account for the thermally
shearing terms. Orthotropic models are easily implemented in an implicit fashion using
an ADI formulation of the finite difference method (e.g., Douglas, 1962) that offers
substantially quicker computations than full-matrix inversions, however, full-matrix
inversions must be performed when modeling any shearing diffusivities using an implicit
finite difference scheme. The shearing terms in the a tensor are not readily implemented
in an ADI formulation, for example, since this method uses operator splitting to divide
each time step into a series of substeps in which only one direction at a time is treated
implicitly (see Press et al., 1992 and Pozrikidis, 1998). It may be easier to implement
finite difference models of the diffusion equation in an explicit formulation when
shearing terms are present in the a tensor. The stability criterion for the explicit problem
is not affected by the shearing diffusivities, as the contribution of these terms negate each
other. The runtime for such a problem would be on the order of twice the explicit
formulation of the orthotropic problem, since modeling of the fully populated 01 tensor
involves computation of three shearing terms on top of the three axial terms already

modeled in the solution of the orthotropic heat diffusion equation.

52

To speed modeling of the thermal response of a randomly oriented anisotropic specimen,
one of the physical axes of the problem, here the out-of-plane axis for a slab specimen,
may be assumed to coincide with a principal thermal direction of the specimen. In such a

case, four of the six off-axis terms would be zero, and the a tensor may be represented as

 

 

0‘11 CL12 0
a: 0121 0I22 0 (3.3)
_ O O (133_

with 002 = (121 due to material symmetry (Powers, 2004). Physically, this tensor would
represent the thermal behavior of an orthotropic specimen with principal diffusivities of
an unknown orientation in the 1-2 plane. Given such behavior, the parameter estimation
routine could be used to seek a solution based on the physical axes imposed here, leading
to the four estimated parameters (11 1. 0122, 0112, and 0133, or, a solution could be sought
based on the principal axes of the specimen, requiring estimation of four different
parameters: the principal thermal diffusivities of the specimen, 01,, 0t", and 01111 (which
here, by definition, is equal to 0133), and B, the angle between the principal and physical
axes in the 1-2 plane. The assumption of one physical axis aligning with a principal
direction is a reasonable one; as for an example, an elastomer rolled into thin plates may
exhibit such behavior, as the manufacturing process for such plates may be shown to
impose anisotropy in the plane of the specimen while leaving the out-of-plane diffusivity
unaffected (Doss and Wright, 2000). As another example, the physical geometry of
lamellar composites favors such an assumption, since the orientation of the lamellae
typically does not favor reorientation of the principal axial diffusivity away from the out-

of-plane axis. Assuming the suggested thermal orientation, the goal would be to design

53

an optimal experiment to resolve the non-zero components in the a tensor without
knowledge of the angle needed to rotate the specimen into the principal orientation and
thus be able to make better measurements of the given a tensor in materials without prior

knowledge of specimen orientation.

Optimization of thermophysical property measurements becomes increasingly
challenging as the complexity of the experiment increases. While 1-D thermal property
experiments such as the flash method (Parker et al., 1961) or guarded hot plate method
(e.g., Bradley, 1943) are simple in design and therefore relatively straightforward to
optimize, experiments used to estimate thermal properties in two or more directions
require well defined heating protocols, boundary conditions and sensor positioning and
thus a more involved optimization analyses. Several thermo-optical methods available to
measure the directional components of the thermal diffusivity tensor include the thermal
wave (Mandelis, 1985), transient grating (Venerus et al., 1999b), and extended flash
methods (Maillet et al., 1990; Doss and Wright, 2000). Each of these methods involves
measuring the response of a specimen to localized heating and each should be designed to
provide the greatest certainty in the property values reported. Design parameters include
the specimen size and shape, the heating area, duration, and magnitude, and the position
and number of the temperature sensors. Prior work has focused on the optimization of
sensor position and heating power and duration for thermal diffusivity measurements on
an orthotropic slab specimen (Davis and Wright, 2005). There, three parameters of
interest (or, 1, 0122, 033) and two nuisance parameters, defined as parameters that must be

estimated along with the parameters of interest in order to solve the inverse problem, but

54

are not valued for themselves (Beck and Arnold, 1977), in this case, incident heat flux
and top face convection, were estimated. Here, the temperature fields are calculated, so
what might be nuisance parameters are assumed known, leaving only the four non-zero
components of the diffusivity tensor to be estimated. The specimen geometry is held
fixed, as are the sensor locations and heating level and duration, which were all based on
results from the prior study. The current investigation focuses on ( 1) an investigation of
the nature of the sensitivity coefficients resulting from a material exhibiting in-plane
anisotropy and random orientation, (2) quantification of the error inherent in assuming
orthotropic behavior for a specimen whose principal thermal axes are rotated an arbitrary
angle from the physical axes, and (3) optimization of the measurement duration given a

specimen of random in-plane anisotropy.

55

3.1 Methods

Estimating thermophysical properties requires repeated solution of the diffusion equation
with the goal of sequentially reducing the difference between the model results and
measured temperatures. For a simulated experiment, data are first generated and
substituted for measurements as input for the inverse problem. Here, an explicit finite
difference model solution is compared to the bottom face temperature history of a
hypothetical specimen as measured at four points. The analytical problem models an
orthotropic planar specimen that has equal in-plane dimensions of length 21 and a
thickness 5 that is much smaller than the in-plane half-width l. The specimen is heated
with a flux go for a finite pulse time over a circular region on the top (2 = 0) surface.
Convection acts over the top and bottom surfaces, while the side surfaces see small
temperature excursions that may be modeled as adiabatic surfaces. The mathematical

formulation of this problem may be stated as

3T =k 32y k BZT 2k 32y +k 32y 3.4
'06 At 11 ax2+ 22 ay2+ 12 axay 33 azz ( )

mOSxSLOSySLOSzS&t>0
Initial Condition: T(x, y, z, 0) = T... (3.5a)

Boundary Conditions:

-k33a%z+h,(T—T,,)=q0 atz=0,forOSx2+y2Sr2,0StStp (3.5b)
h,(T-T,,,,)=k333T8z z=0,forr2Sx2+y2Sl,0StStp (3.50)
h,(T—T,,.,)=k333Taz atz=0,forxSl,ySl,t>tp (35d)

56

8%xz0atx=0andx=l,fort>0 (3.5e)

a%y=0aty=0andy=l,fort>0 (3.51)
h,,(T-T,,)=k3337az atz=8, fort>0 (3.5g)

where k,-,- refers to the ijth component of the thermal conductivity tensor, T... the ambient
temperature, it, the convection heat transfer coefficient of the buoyancy-driven flow at the
top surface, hi, the convection heat transfer coefficient for the stably stratified air layer at
the bottom surface, r the radius of the circular heated region, and t, the heating time.
During parameter estimation, the thermal diffusivity components 0:,- (a;,- = kij lpc) are
estimated instead of the conductivity components kij. Figure 3.2 illustrates the model of
the specimen geometry. The specimen is 40 mm by 40 mm by 3.25 mm thick. The
heated region on the t0p face has a 6 mm radius and is illuminated for 1 s. The
magnitudes of the top and bottom surface convection and heating boundary conditions
are set as known, fixed parameters that are based on values obtained from an
experimental apparatus used to measure the thermal diffusivity of slab specimens.
Temperature sensors are located at four positions on the bottom face of the specimen: one
sensor is at the center of the illuminated area, two of the sensors are located 9 mm from
the center sensor along the coordinate axes defined by the plane of the specimen, and the
final sensor forms the fourth comer of a square 9 mm to a side. The four sensors are
named as shown in Figure 3.2. The finite difference model operates within a parameter
estimation program that is based on a Levenberg-Marquardt iterative nonlinear least
squares fitting routine (Press et al., 1992). This program also calculates show that an 81

x 81 x 41 node mesh produced grid independent temperature histories to within 0.001 K.

57

 

     

corner

y—sensor SGHSOI’
D O

   

 

 

 

 

 

Figure 3.2 Schema of the top surface of the specimen illustrating the illumination area.
Top and bottom faces are subject to specified constant and uniform convection
coefficients. Heating on the top face is of known power and duration; lateral edges are
adiabatic. The temperature sensors are located on the bottom surface of the specimen at
the comers of a 9 mm square that is aligned with the coordinate axes and has one sensor
at the center of the bottom face; the projection of these sensor positions are shown on the
top surface of the specimen. The principal axes of a, x’ and y’, are oriented at an angle B
from the experimental axes; B is assumed unknown at the start of parameter estimation.

Optimizing the duration of the experiment requires a criterion to judge the optimality of
the design. A convenient optimization criterion for inverse heat transfer problems is to
minimize the hypervolume of the confidence interval in the estimated parameters as
measured in the covariance matrix (Beck and Arnold, 1977). Of the several methods

available, D-optimum design is suited for estimation of multiple parameters because the

58

maximization of the D-optimum criterion, the determinant of the Fisher information
matrix A, results in the minimization of the hypervolume of the confidence region about
the estimated parameters (Beck and Arnold, 1977; Mejias et al., 1999; Emery, 2001 ).
The Fisher information matrix is composed of scaled sensitivity coefficients F (F =
p(8T/dp), where p is an estimated parameter and F has units of temperature) that indicate
the relative change in the temperature field at a point with variation in a given estimated
parameter. Comparing different experimental designs requires a non-dimensional

definition of the Fisher information matrix, which may be written as

<0 — 1 .irf[p an“, BEN—1 Yd! (36)
'"'" M4, 3:] ,=0 map", "8p, Tm, '

where M is the number of sensors used (four here), If the total time over which data are

 

 

 

collected, pm and p" the parameters being estimated, T, the temperature at sensor 3 as
calculated by the model using the current estimate for p,, and Tm the highest temperature
reached during the experiment. Note that the first two parenthetical terms inside the
integral are scaled sensitivity coefficients. These terms are already calculated during the
Levenberg-Marquardt parameter estimation, thus calculation of (I), and therefore A,

requires little extra computational time.

To investigate the effects of measurement time and random orientation on the ability to
resolve components of a, temperature histories were calculated for a specimen
possessing an out-of—plane diffusivity of 1.1 x 10'7 m2/s. The axial in-plane diffusivities,
0111 and 022, varied depending on the orientation of the principal axes of a, from a

minimum of or" = 1.317 x 107 m2/s to a maximum of a. = 1.883 x 10'7 mz/s. The

59

corresponding shear diffusivity (002 = 021) varied from zero, where B = 0°, to a maximum
of 0.282 x 10'7 mzls at B = 45°. The magnitudes of the components of a are on the order
of those measured in manufactured PVC (Doss and Wright, 2000) or in loaded
deformable elastomers (LeGall, 2001; Venerus er al., 2001). The values of on 1, (122, and

0112 are related to B and the principal in—plane diffusivities on and an by

(1 +0. (I. ‘11
0L11 =(—I-§—Il]+[I—2£)'COS(B) (3.7)

(1 +0. 0. —0.
0L22 = [I—ZJI-j—[iz—HJ'COSW (3.8)

0. —0. .
a..=[I—21—I)-sm0) (3.9)

Representative values of or; 1. 022, and 0112 are shown in Table 3.1. At the principal
orientation, where B = 0°, the values of the in-plane axial diffusivities are at a maximum
or minimum, respectively, and the shear diffusivity is zero. As B increases, the
magnitudes of the axial diffusivities tend towards the median value of the principal in-
plane diffusivities 011 and or" while the shearing diffusivity increases. When B = 45°, the
shearing diffusivity is at a maximum and the in-plane axial diffusivities are equal to the
median value. For 45° < B < 90°, the shearing diffusivity decreases while the axial
diffusivities again approach maximum or minimum values until B = 90°, where 01 attains
the same numerical values as initially, except that the values of or” and 0122 are reversed.
Temperature histories were generated at a sampling rate of 120 Hz over 150 seconds for
the five values of B shown in Table 3.1. Three thousand data points were used for each

case, resulting in six combinations of sampling rate and measurement duration ranging

60

Table 3.1 Representative values of B, on 1. 022, and 012 used in the analytical study.

 

Test B (1.11 x 107 (mzls) 02 x 107 (mZ/s) (112 x 107 (m2/s)
1 0° 1.883 1.317 0.000
2 15° 1.844 1.356 0.141
3 225° 1.800 1.400 0.200
4 30° 1.741 1.459 0.244
5 45° 1.600 1.600 0.282

 

 

 

 

from 25 s at 120 Hz to 150 s at 20 Hz. Noise was added to each temperature history in
order to simulate the experimental measurements that would normally serve as input for
thermal property estimation. The noise had a Gaussian distribution with zero mean and a
standard deviation of o = 0.01 K, which is on the order of the noise measured using IR
thermography. The resulting thirty sets of data served as input for parameter estimation
where the four non-zero components of thermal diffusivity were estimated, sensitivity
coefficients were calculated for each estimated parameter at each of the four sensor
positions, and A was calculated to provide a relative measure of the amount of

information obtained from each simulated experiment.

Investigation of the error inherent in assuming in-plane orthotropy for samples having a
randomly oriented a required smaller variations in 7 than used in the above tests,
therefore a total of 12 data sets were calculated for 0° S B S 45°. These temperature
histories were sampled at a rate of 30 Hz over a duration of 100 s for a total of 3000 data
points per sensor. These tests were used to develop sensitivity coefficients for each
sensor position as well as to serve as input for parameter estimation. Noise, with the

same attributes as before, was again added to each temperature history. To enforce

61

orthotropy during parameter estimation, the value of the shearing diffusivity was fixed to
zero during the estimation process, leaving only the three axial diffusivities to be
estimated. Once again, A was calculated to serve as a metric of comparison for the

amount of information obtained from each case.

62

3.2 Results

For each orientation of or listed in Table 3.1, the error in the estimated diffusivity, e =
(05- 0to)/010, where as is the estimated diffusivity and org the input diffusivity for the
calculated temperature profiles, was less than 0.1% for the axial diffusivities and less
than 1% for the shearing diffusivity. The fit between the simulated temperature histories
and the inverse solutions generated using the estimated diffusivities was within the
standard deviation of noise in the data. Figures 3.3a and 3.3b show representative
temperature histories and the corresponding inverse solutions for case 1 (B = 0°, where
the principal axes of a. and physical axes of the specimen coincide) and case 5 (B = 45°,
where the shearing diffusivity reaches a maximum), respectively. The inverse solutions
fit the simulated data well as evidenced by residuals (defined as the difference between
the modeled and measured temperature histories) that were computed for each case. As
expected, the residuals oscillate around zero with a magnitude on the order of the noise in
the input data. The x- and y-temperature histories overlap each other since the in-plane
axial diffusivities will be equal for the orientation corresponding to the maximum
shearing diffusivity (Figure 3.3b). The sensitivity coefficients at each of the four sensor
positions in case 1 (Figure 3.4) and case 5 (Figure 3.5) exhibit uncorrelated, nonzero
behavior, with the exceptions of the an coefficient in case 1 and the or, 1 and 022
coefficients at the center sensor in case 5. Recalling the definition of the scaled
sensitivity coefficient F, when the shearing diffusivity has zero magnitude, as it does for
the B = 0° orientation (Figure 3.4), the sensitivity coefficient should be equal to zero over

the entire duration of the experiment. Note how the 012 sensitivity coefficient has a zero

63

 

 
 

x-sensor

center sensor

   

   

 

 

 

 

 

 

     
 

y-sensor
0 I Icomer sensor
_1 . , . , . , f , .
0 20 40 60 80 100
time (s)
5
, ('9)
center sensor
4 -
3 —4
:x: 2 -
®

X-SBDSOI’, y-sensor

HAwAw

   

comer sensor

 

 

 

-1 ' 1 * 1 * 1 ' I ‘
O 20 40 60 80 1 00
time (s)

Figure 3.3 Simulated temperature excess, 9 = T-T... and inverse solutions for (a)

orthotropic test (case 1) and (b) maximum shear test (case 5). Note that for case 5, the x-
and y-sensor temperature histories are identical as a result of 0m and 0122 being equal.

64

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.8 0.16 ......
0.4 I '0'" ...~...“-..... "31:... ' ‘
g ‘ a" " Q 0.08 o.’
V 0.2 v ,9
l— l. c
0 0 0.04 ‘ c;
-0 2 ~ 0 00 .J
-0 4 . . r . -004 . . . .
0 20 40 60 80 100 333 0 20 40 60 80 100
time (s) ......... at 1 time (s)
------- a22
(c) ....... - .12 (d)
6 0.8
0.6 ‘
0.4 " ...........
g ........
,_ 0.2
0.0 emu“.
-o.2 . “New."
-2 . . . . .04 . . . .
0 20 40 60 80 100 0 20 40 60 80 100
time (s) time (s)

Figure 3.4 Scaled sensitivity coefficients for case 1: (a) y- sensor, (b) corner sensor, (c)
center sensor, ((1) x-sensor. The solid line represents F for 033, the dotted line F for or, 1.
the dashed line F for 0122, and the dot-dash line F for (112. Time is not scaled, as the
choice of appropriate scaling parameters is not immediately obvious. For example, if the
specimen thickness and corresponding axial diffusivity (133 are used as scaling
parameters, than 1 second would correspond to a Fourier number of 0.0104; however, if
either of the principal in-plane diffusivities are used with the specimen width as scaling
parameters, 1 second would correspond to a Fourier number of either 8.23 x 10'5 or 1.18
x 10 .

65

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(a) (b)

0.8 0.16

0 6 I'l"..........“~. .. 012 1 '/.—”-~\\

Q4 . ." ........ / \3.‘
2 '0 ‘. ..... A 0.08 .1 '/ ’.'.’ ............. Q:
I: 0.2: 5 // /,.’°’

L— 0.04 ‘ r
0.0 r — . f\
-0.2 ,, 0.00 ‘
-0.4 . . . v .004 r . . .
0 20 40 60 80 100 a33 0 20 40 60 80
time (s) ......... a11 time (s)
------- a22
(C) ------- - 312 (d)
6 0.8
06 4 ....................
4 . ...........
0.4 ............
Q q g .....
,_ 2 L- 0.2
o 0.0 fl...
\.._ 2.- 0.2 “‘~ ....... ,,__ __

-2 . 4 . . . .04 . . . .

0 20 40 60 80 100 0 20 40 60 80
time (8) time (s)

 

 

 

 

 

Figure 3.5 Scaled sensitivity coefficients for case 5: (a) y-sensor, (b) comer sensor, (0)
center sensor, ((1) x-sensor. The solid line represents F for 0133, the dotted line F for on 1.
the dashed line F for 022, and the dot-dash line F for 002.

value at all sensor locations in case 1 (Figure 3.4), while there is a change in the

magnitude of the m2 sensitivity coefficient at each sensor position in case 5, where the

shearing diffusivity is a maximum (Figure 3.5). The amplitude of the 012 sensitivity

coefficient increases at each sensor location as B and the shearing diffusivity increase. In

order to estimate a particular parameter to a high degree of confidence, the magnitude of

the corresponding scaled sensitivity coefficient should be on the order of the temperature

66

 

100

 

100

signal at that point and the behavior of all of the sensitivity coefficients at that location
should be uncorrelated. Over time, large variations in coefficient magnitude are
preferable to relatively constant trends. The magnitude of the largest sensitivity
coefficient at each sensor position is on the order of the temperature signal at that sensor.
The sensitivity coefficients reach their maximum values at different times, which is

preferable when estimating independent parameters.

Assuming thermal orthotropy (i.e., assuming B = 0°, leading to shearing terms of zero
magnitude) for a specimen having principal values of a = [1.10, 1.88, 1.32] x 10'7 m2/s
with randomly oriented in-plane axes results in small errors in the estimated diffusivity
values over the range of in-plane orientation 0° S B S 45° (Table 3.2). Note that the error
in the estimated diffusivity value never exceeds 0.4%. The small error in the axial
diffusivities results from the shearing diffusivity value being at most 18% of the value of
the mean axial diffusivity. This error increases as the magnitude of the shearing
diffusivity increases, e.g., this error increases by an order of magnitude when the
maximum shearing diffusivity is doubled as shown in the first half of Table 3.3.
Increasing the maximum shearing diffusivity to 72% of the median axial diffusivity value
increases the maximum estimated axial diffusivity error to greater than 9% (Table 3.3).
Also included in Table 3.2 is the normalized x2 for each test, which gives an indication of
the quality of fit between the input and inverse model temperature histories. The x2 has

units of temperature squared and is defined here as

N 2
1
X2(Pi)=-A72[Tj “Yj(Pi)] (3.10)
i=1

67

where N is the number of data points, p,- the parameters being estimated, T,- the measured
(simulated) temperatures, and y,-(p,-) the model value corresponding to the jth temperature.
Based on this definition, x2 should have a value on the order of the square of the
magnitude of the noise in the measurements. Assuming that the estimated parameters are
known precisely and that the noise has a Gaussian distribution with a zero mean and a
magnitude on the order of $0.01 K, the value of x2 should be 4 x 10“1 K2. As B increases
for a given specimen, the shearing diffusivity increases and the simulated material
behavior becomes less orthotropic. While the estimates of the axial diffusivities change
little with increasing B, the overall )8 increases, meaning that the quality of the fit
between the temperature histories of the model and the input data decreases. Note in
Tables 3.2 and 3.3 that this error increases with increasing B, corresponding to increasing
shearing diffusivity. This degradation of the overall ability of the model to fit the input
data is manifested in the comer temperature sensor (Figure 3.6). Note how the model
trace no longer lies within the data scatter, indicating that data obtained from the comer
temperature sensor may be used to verify whether the temperature sensor placement on
an anisotropic specimen is coincident with the principal directions of diffusivity, and
thus, whether an orthotropic thermal response may be assumed during the experiment.
These data suggest that when the normalized x2 approximately 20 times the level of the

noise squared (x2 = 84.3 x 104 K2), estimates of the orthotropic components of thermal

diffusivity may be in error by 1% or more.

68

Table 3.2 Error in estimating diffusivity when assuming orthotropy (i.e., neglecting 002)
for an anisotropic material with principal axes oriented by B with respect to the
measurement axes. The maximum shear diffusivity for this material is about 18% of the
median axial diffusivity value.

 

Case (3 e in 6133 (%) e in or” (%) e in (122(%) 22 x104 (K2)
6 0° 0 0 0.08 3.95
7 25° 0 0 0.15 4.24
8 75° 0 0 0.23 6.19
9 11° 0 0 0.22 8.65
10 15° 0 -0.05 0.22 12.5
11 18° 0 -01 1 0.22 16.0
12 225° 0 -0.17 0.14 21.7
13 26° 0 -0.17 0.07 26.4
14 30° -0.09 -0.23 -007 31.9
15 35° -0.09 -035 .007 37.0
16 40° -0.09 -03 -019 41.1
17 45° -0.09 -0.38 -025 42.1

 

 

 

 

 

Table 3.3 Error in estimating diffusivity when assuming orthotropy for an anisotropic
material with principal axes oriented by B with respect to the measurement axes. The

maximum shear diffusivity for this material is about 35% of the median axial diffusivity
value for tests 18-22 and about 64% of the median axial diffusivity value for tests 23-27.

 

Case (3 e in (133 (%) e in or“ (%) e in 6122 (%) x2 x104 (K2)
18 0° 0.00 0.00 0.10 3.98
19 15° .0.09 -057 0.00 40.3
20 225° -027 -1.10 -075 84.3
21 30° -0.55 -1.65 -152 134.
22 45° -0.73 -231 -225 191.
23 0° 0.00 0.00 0.00 4.02
24 15° -0.64 -1.49 -307 127.
25 25° -l.64 -359 -7.52 416.
26 32° -2.45 -5.62 -9. 19 670.
27 45° -3.18 -8.88 -8.81 924.

 

 

 

69

 

 

Optimizing the duration of the experiment for a given (1 having principal values of [1.10,
1.88, 1.32] x 10'7 mzls and a heating time of l 8 indicates that the maximum A is obtained
at 75 8, assuming that the amount of sampled data is held constant between tests. Figure
3.7 shows the variation of normalized A as a function of measurement duration. Note
that A/Amam is maximized for an measurement duration of between 75 and 100 s. The
only case with a distinct maximum in A/Amax is for B = 0°, where the maximum lies at
100 s. Non-normalized values of A led to a maximum at 75 s for the orientation angle [3

= 45°.

70

 

 

    

 

      
      

 

 

 

 

1.0
x-sensor, y-sensor
-' «f ,J'xa.§£t...4.w.u.s {’19.
0.8 7 1
0.6 ‘
A
it
Q corner sensor 1 _
0.2 7
0.0 ‘ r’" ‘ ,
'O.2 ' l ' I r l ' 1
0 20 40 60 80 100

time (s)

Figure 3.6 Input data and model tits at the x-, y-, and comer sensors for case 5 (B = 45°)
where orthotropic behavior is assumed during parameter estimation. The solid lines
represent the model fits for the simulated temperature history at the indicated sensor. The
poor fit between the data and the model temperature histories at the comer sensor is due
to the assumption of no shearing diffusivity. The temperature responses at the x- and y-
sensor locations overlap since the x- and y—diffusivities are equal.

71

 

 

 

 

1.0 7 ' 9
U
0.8 ‘
l
0.6 7
<8
i 0
a 0.4 7
5 O O
0.2 "
Q?
0.0 “ ' O
T l ' I ' I T l 7 l r
0 25 50 75 100 125 150

time [s]

Figure 3.7 Normalized A as a function of measurement duration for cases 1 through 5.
Empty circles represent data taken for B = 0°, triangles B = 15°, squares B = 22.5°,
diamonds B = 30°, and filled circles B = 45°.

72

3.3 Discussion

For the cases where all four non-zero components of a were estimated, the maximum
error in the estimated values of the axial components of a was within 0.1% of the input
value, while the maximum error in the estimated shear diffusivity was within 1% of the
input value. The error in the estimated axial diffusivities is on the order of that seen by
Graham et al. (1999a), where the estimated parameters were the principal diffusivities
and the orientation angle B. In that study, the maximum error in the estimated
diffusivities was 0.2% of the input value. Here, the shear diffusivity error increased from
0.1% to nearly 1% for very small B, where the value of the shearing diffusivity
approached zero, magnifying the effects of small estimation errors. Inspection of F at
each sensor location for the four estimated components of (1 showed uncorrelated
behavior for each non-zero component of (1, provided that data from all four locations are
taken into account. Neglecting the behavior of the om and an sensitivity coefficients at
the center sensor, each sensitivity coefficient reaches its maximum at a different sensor
position, and it is at that position where the most information is obtained concerning the
parameter associated with each coefficient. For example, one expects that much of the
information about 0133 is gained from the center sensor during the first 20 to 30 seconds of
the experiment (Figure 3.4). Measurements on PVC conducted via a three-sensor format
have shown this to be the case. The behavior of the sensitivity coefficients shows that the
x- and y-sensor locations provide most of the information used to determine an and (122.
Because the or“ and org; sensitivity coefficients have nearly the same behavior at the

center sensor position (or even identical behavior when their respective diffusivities are

73

equal, as shown in Figure 3.5a), it would be difficult to determine each parameter
independently using that sensor position alone, even though the magnitudes of these
coefficients are higher at the center position than any of the other sensors. The responses
of these sensitivity coefficients are clearly independent, however, at either the x- or y-
sensor locations (Figures 3.5c and 3.5d). This result indicates the importance of using the
appropriate sensor locations to estimate multiple parameters in thermophysical property
measurements. The comer sensor location delivers the most information about 0112,
because the (112 sensitivity coefficient has its highest value there. Additionally, the
magnitude of the 0112 sensitivity coefficient is on the order of the other sensitivity
coefficients at that location only. Such behavior indicates that the four-sensor array
described in this study provides adequate information for resolving the four components
of 0. that are of interest here. Furthermore, the location of the sensors affects the ability
to estimate the components of a. If the two in-plane axial diffusivities are to be the only
estimated parameters, and the center (Figure 3.5a) and corner (Figure 3.5b) sensor
locations were used to estimate these parameters, the correlated behavior of the 0m and
an sensitivity coefficients at both locations would make the independent estimation of
each parameter difficult, or even impossible. Because specimen orientation may be
unknown prior to the experiment, sensor placement on orthogonal axes may ensure the

best chance of obtaining uncorrelated sensitivity coefficients for the parameters of

interest.

Estimating uwhen orientation along principal axes has been assumed for a specimen that

is randomly oriented resulted in errors in the estimated axial diffusivities below 1%

74

irrespective of the specimen orientation when the maximum shearing diffusivity of the
specimen was about 18% of the median axial diffusivity. This error increases as the
magnitude of the shearing diffusivity increases, as shown in Table 3.3, however, the
magnitude of anisotropy typically seen in experiments on elastomers is closer to that of
the low-shearin g diffusivity cases listed in Table 3.2. The in-plane to out-of-plane
anisotropy ratio calculated in residually stressed PVC was 1.25 (Doss and Wright, 2000),
while recent work on loaded elastomers reported anisotropy ratios below 1.2 for stretches
of up to twice the original specimen length (Venerus et al., 2001). Here, the ratio
between principal in-plane diffusivities was 1.43 for the low shear cases (Table 3.2). The
anisotropy ratio will increase as the maximum shearing diffusivity increases, provided
that the median axial diffusivity is held constant (Table 3.3). Thus, while a high level of
anisotropy may affect the accuracy of results obtained under assumptions of orthotropy,
such levels are unlikely to be encountered except in composites that are manufactured to
have such properties, therefore, orthotropy may be assumed during property estimation
for materials possessing low levels of anisotropy with little loss of accuracy. Still, if high
levels of accuracy are needed during property estimation, it may be advantageous to
measure 0. at orientations that are not aligned with the principal axes of or once a
specimen is determined to exhibit in-plane anisotropy. Figure 3.8 shows the variation of
Awith B for a given a. Note that A increases as B increases, despite the number of data
points being the same for each B. While the error in the estimated diffusivities is always
less than 1% for the lowest shearing diffusivity cases (Table 3.2), suggesting that
specimen orientation had little effect on parameter estimation, the amount of information

obtained during each test as quantified by A increased by at least an order of magnitude

75

when the specimen orientation was rotated away from the principal thermal axes.
Because maximizing A minimizes the hypervolume of the confidence region of the
estimated parameters (see Mejias et al., 1999), A is inversely related to the uncertainty in
the estimated parameters, therefore, parameter estimation experiments performed on
anisotropic specimens oriented at B = 45° with respect to the principal axes of a should
yield more precise results, if the orientation is known prior to testing. Given that the
boundary conditions here were assumed perfectly known, the effects of additional
information obtained at the off-axis orientations may help increase accuracy in actual
experiments where only approximations of boundary conditions may be known or

significant uncertainties exist in geometric parameters such as sensor location or

specimen thickness.

Investigation of the optimal measurement duration for a particular or indicated that the
optimal duration was between 75 s and 100 s for each orientation angle B that was
studied. A prior study on orthotropic materials possessing identical in-plane axial
diffusivities showed that the optimal measurement duration occurred when the in-plane
temperature histories reached their maximum values (Davis and Wright, 2005). Other
prior studies have argued on physical grounds that the optimal measurement duration
coincides with the maximum in-plane temperature rise. Doss (1998), for example, found
that the estimated diffusivities om of A181 304 stainless steel varied by more than 5%
when the experiment duration was nearly half the optimal value, but the estimated values
were within 2% 0f the reference value when the optimal measurement duration was used.

l/eGall (2001) found that the estimated properties of plate glass varied by less than 0.5%

76

when the optimal measurement duration was used, but when the measurement duration
was altered by i50%, the scatter in the estimated diffusivities increased to 5%. The x-
and y-sensor temperature histories in this particular case reached their maximum values
between 75 s and 100 s for the B = 0° orientation. Since the alignment of the specimen
axes with the principal axes of (1 leads to orthotropic behavior at this orientation, these
results resemble those from prior studies, but with the extension that the orientation angle
has little effect on the optimal measurement duration. Based on the results of the prior
study, one might expect the optimal measurement duration to be the same as the longest
sensor rise time, which would occur in this study at the comer sensor, but actually the
optimal measurement duration occurs at an earlier time. The interaction of the
components of the diffusivity tensor with the fixed boundary conditions renders the

problem complicated enough that computational simulation is required to find the

optimal experimental time.

77

 

10'3 i
10'4 ‘3 o
10'5
10'6
10'7
10'8

10'9 “a

 

 

1 0‘1 0 1 ' l . 1 ' 1 ' I '
O 10 20 30 4O 50

Figure 3.8 Variation of A for an or having principal values of [1.10, 1.88, 1.32] x 10'7
m /s. The amount of information obtained from a measurement increases as the

specimen orientation shifts away from the principal axes of a, Note that for B = 0° and
25°, A is six and two orders of magnitude smaller, respectively, than A at B = 45°.

78

3.4 Closure

An analytical study was performed where the four non-zero components of the diffusivity
tensor (1 were estimated for a material possessing random in-plane anisotropy on the
order of certain manufactured or mechanically loaded elastomers. Sensitivity
coefficients were calculated for each of the non-zero components of the a tensor at four
sensor locations. The amplitude of the sensitivity coefficients suggest that the signal
received at each sensor is sufficient for the estimation of the parameters, while the
character of the coefficients over time suggests that the parameters are independent.
Small values of in-plane anisotropy lead to negligible (< 0.4%) errors in the estimated
parameters when the shearing diffusivity is neglected. Larger differences in anisotropy
require the estimation of the shearing diffusivity, which is possible with the addition of a
fourth sensor to the extended flash diffusivity method (Doss and Wright, 2000). This

sensor may be used to verify assumed thermally orthotropic behavior in materials.

79

APPENDICES

80

APPENDIX A

NOTES CONCERNING IR THERMOGRAPHY

81

The use of IR thermography as a temperature measurement technique for extended flash
thermal diffusivity measurements is particularly desirable due to the flexibility it provides
the experimenter for positioning of the temperature sensors. Additionally, thermo-optical
methods of temperature measurement make no contact with the specimen surface that
could disturb the resultant temperature field, thus eliminating a source of error in the
experiment. However, the use of IR thermography presents its own challenges that must
be considered if such a method is to be used in the lab. Presented here are suggestions
concerning (a) the general use of IR thermography for temperature measurement and (b)
specifically, the use of our own IR camera in the data acquisition and subsequent

parameter estimation procedure.

The camera used for making experimental temperature measurements in these studies is
an Indigo Merlin Mid-IR (3 — 5 pm) model that is Stirling-cycle cooled. The camera uses
a 30 um x 30 um indium antimonide (InSb) detector to record a 320 x 256-pixel field at a
frame rate of up to 60 Hz. The camera lens currently in use has a 50 mm focal length, a
field of view of 11° x 8°, and an instantaneous field of view of 0.6 milliradians. The
response time of the detector, as measured in cameras of similar wavelengths, is on the
order of several microseconds (Torres 42! al., 1990). The camera comes factory calibrated
to measure temperature over the range of 0-350 °C, but can be adjusted via software to
handle temperatures up to 2000 °C if necessary. Typically, the temperature range of
interest for our measurements is 20-40 °C. The noise seen in spot camera measurements
is approximately $0.01 K, about an order of magnitude less than that seen when using E-

type thermocouples with commercial calibration. The camera pixel size depends on the

82

Optics as well as the distance between the specimen and lens and is 210 um square in the
current test configuration. The camera captures full-field data at a constant frame rate of
60 Hz. Using the accompanying software, it is possible to capture data at slower rates by
saving every 11th frame, leading to successive capture rates of 30 Hz, 15 Hz, 12 Hz, 10 Hz,
etc. Slower capture rates are preferable for lower diffusivity materials such as PVC,
polyurethane rubber, or plate glass. The maximum rate of 60 Hz is sufficient for
capturing data from a 1.5 mm thick AISI 304 stainless steel plate, but 3-D flash tests on
such a specimen yield only about 350 data points if the data is sampled until the off-axis
temperature sensors (assumed to lie ~8-10 mm from the center sensor) reach their
maximum temperature (~6 s @ 60 Hz). The center-sensor temperature rise occurs in 2-3
s for materials with diffusivities on the order of steel; should the desire to examine a
material with a higher diffusivity arise, a camera with a faster frame rate should be used
for data capture. As noted in this dissertation, there is a random offset between the laser
and camera triggers due to the fact that the camera does not synchronize precisely with
the laser trigger. For elastomers, this offset is negligible, causing little to no change in
any estimated parameters, but for higher diffusivity materials like steel, several tests must
be run and subsequently averaged to arrive at an appropriate value for the out-of—plane
diffusivity. As for sensor positioning, the act of moving the off-axis sensors further away
from the center point results in a slower rise time, allowing more data points over the
temperature history, but also a lower maximum signal and a poorer signal to noise ratio.
It is possible that the constant convection boundary condition modeled on the top surface
(and estimated as a nuisance parameter) breaks down near the edges where a convective

plume would likely pull cool air in from the surroundings. As a rule of thumb, the outer

83

5-7 mm of the specimen probably shouldn’t be used for data sampling with the finite
difference models included in this dissertation. Additionally, the specimen holder itself
may cause parasitic heat losses. Originally, one of the sensor positions used in the study
was located between the illuminated area and the specimen holder (which clamps to the
specimen at one point via a screw). The metal holder was found to affect the temperature
response at this sensor location, causing a lower temperature signal than what was
expected. The responses at the other off-axis locations yielded an isotropic response. It
is recommended that sensor positions always be located away from the point where rigid

specimens are clamped to the device.

The detector records the infrared luminosity of the target and converts that number to a
temperature value. Unfortunately, the target luminosity can be biased by reflected
infrared light from elsewhere. The entire testing apparatus has been covered with a tarp
that is secured to a rigid frame to minimize intrusion of light from outside sources, but
light still leaks in from around power cords, and the camera itself generates heat which
can unevenly warm the apparatus and surroundings. Surface reflectivity effects were
seen in the full field images when testing barren aluminum and steel plates, and minor
effects were also seen when testing PVC plates. To counteract the reflectivity, flat black
stove spray paint was used to coat the top (illuminated) and bottom (interrogated)
surfaces of the specimens. This not only reduced the surface reflectivity of the specimen,
but also the transmissivity of the laser light through translucent specimens such as glass.
The thickness of the paint is on the order of microns, several orders of magnitude smaller

than the specimen thickness, and thus has a negligible effect on the estimated parameters

84

of the experiment. In prior experiments, a similar paint was used to eliminate
transmission of laser light through a glass specimen while using a fixed thermocouple
probe for temperature measurement. At the time, it was believed that the laser energy
was completely absorbed at the coating surface, but IR camera tests on the plate glass
specimen showed that the laser light penetrates through both painted surfaces resulting in
a small but noticeable signature in the center sensor temperature history that ceases when
the laser pulse ends. Thus, the paint is opaque to visible light, but lets infrared light
through. Multiple coats of the paint do not remedy the problem. A similar signature is
not seen in the center sensor temperature history of any metal or elastomeric specimens.
Should the desire to test a translucent or transparent material arise, another material,
preferably one that may be applied via spraying so as to minimize the amount of
deposited material, must be used. Since the stove paint appears opaque but actually
transmits a fraction of the incident infrared light, possible coatings for the specimen
should be tested using the camera to check for a signature from the laser. The laser
should be tuned to its highest power for a 1-2 s length pulse time, depending on the
thickness and estimated properties of the material. If there is transmission, a single
sawtooth-shaped spike will appear in the center temperature sensor temperature history at
the end of the flash time. The spike seen during glass measurements was roughly one to
two orders of magnitude greater than the ambient noise. If transmission is absent, the
temperature rise at the center sensor should be smooth (within the bounds of the

experimental noise).

85

When data is collected, it is saved in a Dynamite movie file. This file can then be saved
frame by frame. These single frames may be opened in the Thermogram software, and it
is from here that temperature histories at specific locations may be pulled using the Plot
tools. After the appropriate sensor areas are selected, temperature histories may be pulled
from each frame and displayed on a visual plot. This data may then be saved as an
ASCII file, and this file is used as input for the parameter estimation program. Currently,
the Marquardt program reads in 3-sensor ASCII data. The format of the data line from
the ASCII file is “date, time, sensorl, sensor2, sensor3”, however, the file formatting that
the camera uses depends on the time at which the camera is being used -— different times
of day lead to different spacings between the time and the temperature values, and thus
when running news tests on a specimen, the data files that are read in must be checked to
ensure that the correct formatting is used for the input statements. The parameter
estimation program pauses after reading in the data file, giving the experimenter a chance
to check the input data, which is echoed to the “dataout” file. If the data file is not read
properly, negative temperatures will appear, or the temperature response at a given sensor
location will suddenly jump up and down by tens of degrees. To correct this, the
FORMAT statement that is used to read in the data inside the “readit” subroutine in the
Marquardt program must be altered. Comment lines inside the program list the changes

necessary to correct for several common formatting errors.
To approximate point-wise temperature measurement, the Thermogram software is used

to record data at three points using the “rectangle” tool to record average temperatures

over a 3 X 3-pixel area. This sensor area is of sufficient size that the detector will capture

86

the true temperature of the position of interest yet minimize noise. This 3 x 3-pixel array
produces a sensor area of about 0.6 mm by 0.6 mm. While this area is smaller than the
surface area of the E-type thermocouples previously used for the fixed-position
temperature probe, there is still some uncertainty in the position of the sensor. Sensor
positions should be determined using a reference material with known thermal properties.
Setting the thermal diffusivities as known quantities, estimate the off-axis sensor
positions, incident heat flux and top surface convection coefficient using the Marquardt
program. Once the sensor positions are confirmed, materials of unknown thermal
properties may be tested. The clamp inside the testing device is set such that the bottom
surface of the (rigid) specimen is always the same distance from the camera lens; as long
as the camera and clamp are not moved, the pixel dimensions will remain the same and

sensor positions do not need to be recalculated.

87

APPENDIX B

ALTERNATING DIRECTION IMPLICIT (ADI) FINITE DIFFERENCE

PROGRAM FOR MODELING ORTHOTROPIC HEAT CONDUCTION IN A

SLAB SPECIMEN

88

The following pages contain FORTRAN code for an alternating direction implicit (ADI)

finite difference program that models the diffusion equation and assumes thermal

orthotropy in a slab specimen. The thermal diffusivity tensor (1 that allows an orthotropic

response may be expressed as

 

111] 0 O _
9= 0 922 0 (8.1)
_ O 0 (133a

 

The simulated specimen boundary conditions consist of: insulated edges, convection on
the top and bottom surfaces (of different magnitude), and a top surface heat pulse
approximately 5.3 mm by 5.3 mm with constant magnitude and finite duration. The

estimated parameters include the three axial components of a, a scaled heat flux

coefficient, and a scaled convection coefficient for the top surface. The bottom surface
convection coefficient is assumed known. Most of the experimental and geometric
parameters are adjusted via the setup subroutine that precedes the actual model routine.

The output consists of an array containing the temperature history of the simulated

specimen at three points on the bottom surface.

89

O ()(5 O ()(5 O ()(5 0 (3(3 0 C)

O

Implicit ADI model for IR camera data reduction.
Model is used with Marquardt program to solve IHCP
for a specimen of unknown thermophysical properties.

This monolayer model has a rectangular, uniform surface heat

flux on the top surface of specimen. All side boundaries are

insulated. The top and bottom surfaces have convection BCs implemented
due to unstably stratified air heating on the top and stabily stratified

heating on the bottom (Htop > Hbottom). The bottom convection coefficient

is assumed known.
*******************

this routine sets up fixed parameters for the problem to be solved.
I wanted to set all the variable in one program as opposed to looking
all over creation to change the problem of interest.

subroutine setup(IMAX,STOPITER,XAPT,YAPT)

CHANGE model size here
NMAX should equal NX times NY times NZ for this prog.
useful combinations: (21,21 ,21,9261) (41,41,41,68921) (61 ,61 ,41 , 152561 )

integer NMAX,NX,NY,NZ
parameter (NX = 41)
parameter (NY = 41)
parameter (NZ = 41)
parameter (NMAX = 68921)

real(8) zlength,xlength,ylength,xratio,yratio,ft

real(8) dz,dx,dy,QIB,QJB,fliz,ALPHAZ,ALPHAX,ALPHAY,QEST,DHZ
real(8) pcspec,flashtime,freq,dt,BOTCONV

real(8) endtime,timecalc,STOPITER,scalz,scalx,scaly,SURFABS

real(8) tdiff(7),tdiffo(7),dx2i,dy2i,d22i

real(8) XAPT,YAPT,lxtc,lytc,xcorr,ycorr,Tmax

integer QI,QJ ,nreada,dataswitch,PN X,PNY,PN Z,NTOTAL,N POINTS
integer FLASH,IMAX,x(5),y(5),z(5),nread,XTC1,YTC1,XTC2,YTC2
integer noiseswitch,sensors

character*12 fname

common blocks

common /GUESS l/ tdiff,tdiffo

common IOTPT/ x,y,z,XTC1,XTC2,YTC1,YTC2,xcorr,ycorr

common lPOSlT/ PNX,PNY,PNZ
common /RFU QI,QJ ,QIB,QJ B

90

0000000

0000000

0

common IREADITT/ nreada,dataswitch,fname,scalz,scalx,scaly,fhz
common IGEOM/ dx,dy,dz

common IGEOM2/ dx2i,dy2i,d22i,dt

common /ETC/ NPOINTS,NTOTAL,FLASH,endtime,pcspec
common ITIME/ timecalc

common IFOO/ nread

common lFLASHIT/ flashtime

common /ZL/ zlength

common IOUT2/ noiseswitch

common IFISI-I/ sensors,Tmax

*model specimen property parameters (initial guesses)*

3-D thermal diffusivity, heat flux, convection
ALPHAZ = 1.25D-07

ALPHAX = 1.25D-07

ALPHAY = 1.25D-07

QEST = 0.031

QEST = 0.02

DHZ = 1.25D-05

ALPHAZ = 4.1D-06
ALPHAX = 4.1D-06
ALPHAY = 4.1D-06

QEST = 0.01

DHZ = 1.0D-05

DHZ = 0.0

assume that pc of the specimen can be found thru DSC.
glass

pcspec = 2.268D+06

pvc

pcspec = 1.540D+06

steel

pcspec = 3.866D+06

rear face convection -- assumed known

to reset to insulated boundary, zero out BOTCONV
BOTCONV = 5.0

BOTCONV = 0.0

*model run parameters*
sampling frequency, # of data points, overall calc. time

fhz=15.
fl'lZ = 30.

91

fhz = 60.

nreada = 1190
timecalc = 18.01
flash time
flashtime = 1.0
flashtime = 10.0
flashtime = 1.31

flashtime = 2.2
dataswitch: 1 = IR camera (3-spot) data, 2 = TC data, 3 = contrived test data

4 = forward problem
dataswitch = 1
noiseswitch: 0 = no noise, 1 = add noise
this should be set to "l" for the forward problem only!
noiseswitch = 0
input (data) file name.
fname = "st04c.asc"

*location of 'output' nodes (in mm)*

lxtc = 20.0
lytc = lxtc
lxtc = 17.26
lytc = 16.23

*nodal geometry parameters*

input dimensions in millimeters

the model is nondimensionalized by the specimen thickness
zlength = 3.355

zlength = 2.25

zlength = 1.50

xlength = 20.
ylength = 20.
zlength = l.
xlength = 1.
ylength = 1.
sensors = 3

scaling parameters for TC data

scalz = 0.00

scalx = 0.00

scaly = 0.00

surface absorption only for this model (for now)
SURFABS=1.00

92

0

iteration parameters (max. # of iterations; convergence criterion for alphaz)

IMAX=50
STOPITER= 1 .0D- 10

nondimensionalization

xratio = xlength/zlength
yratio = ylength/zlength

dz = 1./dfloat(NZ- 1)

dx = xratio/dfloat(NX- 1)
dy = yratio/dfloat(NY- 1)
dx21 = l./(dx*dx)

dy21= l./(dy*dy)

dzZi = 1./(dz*dz)

write(*,*) xratio, yratio, dz, dx, dy

heat flux/aperture parameters. The QI/QJ terms specify the highest node
to which the unit heat flux will be applied.

XAPT = 5.85

XAPT = 5.57

YAPT = XAPT

QI = nint(XAPT/(dx*zlength))
Q] = nint(YAP’F/(dy*zlength))
QI = NX

QJ = NY
these terms add a portion of the total flux to the boundary of

the heated surface -- needed if the actual geometry doesn't line

up with the descretized geometry.
QIB = (XAPT-(float(QI)*dx*zlength)+(0.5*dx*zlength))/(dx*zlength)

QJB = (YAPT-(float(QJ)*dy*zlength)+(0.5*dy*zlength))/(dy*zlength)

QIB = 0.
QIB = O.
if((QI.gt.NX).or.(QJ.gt.NY)) pause 'flux defined area too big'

Output nodes are, by definition, usually on the bottom face of the specimen.
Currently, (1) = top face center, (2) = bottom face center,
(3) = xTC, (4) = yTC, (5) = center bottom outside edge of specimen.

points x(3) and y(4) will be written over for the inter-nodal back face points

x(1)=l
x(2)=1

93

CO

x(3) = O
x(4) = 1
x(5) = NX

y(1)= 1
W) = 1
y(3) = l
y(4) = 0
W) = 1

2(1) = 1

2(2) = NZ
z(3) = NZ
2(4) = NZ
2(5) = NZ

PNX = NX
PNY = NY
PNZ = NZ

NTOTAL = PNX*PNY*PNZ

if (NT OTAL .ne. NMAX) pause 'NTOTAL does not equal NMAX'
if (int(timecalc*fhz).gt.nreada) pause 'check nreada -- too small'

freq = fhz*flashtime

NPOINTS = int(timecalc*fl12)

dt = ( l .lfreq)
endtime = dfloat(dt*NPOINTS)*flashtime
FLASH = int(freq)

ft = float(FLASH)/fhz

this section sets the estimated parameters into the array that
will pass back and forth between the model and the MO reduction

scheme.

tdiffl I )=ALPHAZ
tdiff(2)=ALPHAx
tdiff(3)=ALPHAY
tdiff(4)=QEST
tdiff(5)=DHZ

tdi ff(6)=SURFABS
tdiff(7)=B OTCONV

94

this section takes care of inter-nodal IR camera query points

lxtc = lxtc/zlength
lytc = lytc/zlength

XTCl = int(lxtc/dx) + l
YTCl = int(lytc/dy) + 1
XTC2 = XTCl + 1
YTC2 = YTCl + 1

xcorr = l. - ((lxtc - (XTC1-1)*dx)/dx)
ycorr = 1. — ((lytc - (YTC1-1)*dy)/dy)

Tmax = 0.

return
end

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

C

this is the body of the bilayer model program

subroutine model

to alter the size of the model, update the parameters NX,
NY, NZ and NMAX. To facilitate this, search through the
program for the word "CHANGE" (as in model size).

CHANGE model size here
NMAX should equal NX times NY times NZ for this prog.

integer NMAX,NX,NY,NZ
parameter (NX = 41)
parameter (NY = 41)

parameter (NZ = 41)
parameter (NMAX = 68921)

*variable declarations*

integer N TOTAL
integer STEP,COUNT,FLASH,NPOINTS

integer position,i,j,k,noiseswitch,idum

real(8) flashtime,fluxswitch,endtime,kspec

real(8) dt,dx,dy,dz,dx2i,dy2i,dz2i,time
real(8) Fo,G 1 ,G2,Biot,q,Fol ,F02,Fo3,Bind,qnd,nBind

real(8) AX,AY,AZ,BX,BY,BZ,CX,CY,CZ
real(8) zlength,pcspec,backbiot,bBind,anind

95

O

real(8) T(NX,NY,NZ),Tout(50000,5),noise

real(8) al(NMAX),b1(NMAX),cl(NMAX),r(NMAX)
real(8) a2(NMAX),b2(NMAX),c2(NMAX),u(NMAX)
real(8) a3(NMAX),b3(NMAX),c3(NMAX)

real(8) tone(NX,NY,NZ),ttwo(NX,NY,NZ)

real(8) RA(9),RB(9),RC(9),RD(9),RE(9),RF(9),RG(9)
real(8) RH(9),RI(9),RJ(9),RK(9),RL(9),RM(9),RN(9)
real(8) RO(9),RP(9),RQ(9),RR(9)

real(8) atry(7),tmod(50000)

*** common blocks ***

NOTE: temperature arrays are scaled, dimensionless
temperatures. Do not treat T as T(deg. C)!
NOTE: lengths are scaled too. Time is scaled!

common blocks

common [TEMPERATURE/ T,tone,ttwo,Tout
common /PROPS/ Fo,G1,G2,Biot,q,qnd,Bind,nBind

common IRVECT/ RA,RB,RC,RD,RE,RF,RG,RH,RI,RJ,RK,RL,RM,RN,RO
common IRVECT2/ RP,RQ,RR,bBind,anind

common IFLASHIT/ flashtime

common /GEOM/ dx,dy,dz

common /GEOM2/ dx2i,dy21,d22i,dt
common /ETC/ NPOINTS,NTOTAL,FLASH,endtime,pcspec

common lMAINPROG/ tmod,atry

common /ZL/ zlength
common /OUT2/ noiseswitch

non-dimensional (specimen) parameter construction

Fo = (atry(l)*flashtime)/((zlength/1000.)*(zlength/1000.))

G1 = atry(2)/atry(1)

G2 = atry(3)/atry(l)

G1 = 1.

G2 = l.

kspec = atry(l)*pcspec

Biot = (atry(5)*(zlength/1000.)*pcspec)/kspec

q = (atry(4)*pCSpec)
backbiot = (atry(7)*(zlength/l000.))/kspec

*** body of program ***
from here to the "THIS IS THE LOOP" comment is only called once.

96

C

C

*** construct known coefficients ***

initialize "building blocks" for specimen

F01 = F0*Gl *dt*dx21
F62 = Fo*G2*dt*dy2i
F63 = Fo*dt*d22i

Bind = (2.*Fo*Biot*dt)/dz
nBind = -l .*Bind

bBind = (2.*Fo*backbiot*dt)/dz
anind = -l.*bBind

qnd = (4.*dt)/dz

LHS coefficients
internal nodes -- specimen

AX = -1.*Fol
AY = -l.*Fo2
AZ = -1.*Fo3

BX = 2.*Fol + 2.
BY = 2.*F02 + 2.
B2 = 2.*Fo3 + 2.

CX=AX
CY=AY
CZ=AZ

below are RHS coefficients
first, specimen coefficients

RB(1) = Fol
RB(2) = Fol
RB(3) = Fol
RB(4) = 2.*Fol
RB(5) = 2.*Fol
RB(6) = Fol
RB(7) = 2.*Fol
RB(8) = Fol
RB(9) = Fol
RC(l) = 2.*F02
RC(2) = F02
RC(3) = F02

97

RC(4) = F02
RC(S) = F02
RC(6) = F02
RC(7) = 2.*F02
RC(8) = 2.*F02
RC(9) = F02

RD(l ) = 2.*Fo3
RD(2) = 2.*Fo3
RD(3) = F03
RD(4) = 2.*Fo3
RD(5) = F03
RD(6) = F03
RD(7) = F03
RD(8) = F03
RD(9) = F03

do i = 1,9

RA(i) = 2. - 2.*RB(i) - 2.*RC(i) - 2.*RD(i)
RE(i) = 0.
RF(i) = 0.
RG(i) = 0.
RH(i) = 0.
RI(i) = 0.
R] (i) = 0.
RK(i) = 0.
RL(i) = 0.
RM(i) = nBind
RN(i) = 0.
RO(i) = 0.
RP(i) = anind
RQ(i) = 0.
RR(i) = 0.

enddo

RE(2) = -2.*RB(2)
RB(3) = -2.*RB(3)
RB(5) = -2.*RC(5)
RB(6) = -2.*RC(6)
RE(8) = -2.*RD(8)
RB(9) = -2.*RD(9)

RF(2) = RB(2)
RF(3) = RB(3)
RG(5) = RC(S)
RG(6) = RC(6)

98

0

RH(8) = RD(8)
RH(9) = RD(9)

RI(3) = -2.*RC(3)
RI(6) = -2.*RD(6)
121(9) = -2.*RB(9)

RK(3) = RC(3)
RL(6) = RD(6)
RJ (9) = RB(9)

convection coeffs.

RM(1)= 2.*nBind
RM(2) = 2.*nBind
RM(4) = 2.*nBind

RN(8) = nBind
RN(9) = nBind

RO(6) = nBind

RP(l) = 2.*anind
RP(2) = 2.*anind
RP(4) = 2.*anind

RQ(8) = anind
RQ(9) = anind

RR(6) = anind

construct coefficient vectors -- starting with internal nodes
note that the RHS current node is constructed here as well for
the internal current node ONLY.

doi=1,NMAX
a1(i)=AX
b1(i)=BX
cl(i)=CX

a2(i) = AY
b2(i) = BY
c2(i) = CY

a3(i) = AZ
b3(i) = 82

99

c3(i) = CZ
enddo

construct A and C vectors on the boundaries

top comers

a1 (position(1,1,1 )) = 0.
cl (position(],l, 1 )) = 2.*CX

a1 (position(NX,], 1)) = 2.*AX
cl(position(NX,1,1)) = 0.

a1 (position(1,NY, 1)) = 0.
cl(position(1,NY,1)) = 2.*CX

a1 (position(NX,NY,1)) = 2.*AX
c 1 (position(NX,N Y, 1 )) = 0.

a2(position(1,l ,1)) = 0.
c2(position(1,l,1)) = 2.*CY

a2(position(NX,] ,1)) = 0.
c2(position(NX,],1)) = 2.*CY

a2(position(l,NY,1)) = 2.*AY
02(position(1 ,N Y, 1 )) = 0.

a2(position(NX,NY,l)) = 2.*AY
c2(position(NX,NY,1)) = 0.

a3(position(] ,1 ,1)) = 0.
c3(position(l,1,1)) = 2.*CZ

a3(position(NX,1,1)) = 0.
c3(position(N X, 1 , l )) = 2.*CZ

a3(position(] ,NY,1)) = 0.
c3(position(1,NY,l)) = 2.*CZ

a3(position(NX,NY,l)) = 0.
c3(position(NX,NY,l)) = 2.*CZ

100

c top and bottom edges

doj=2,NY-l

enddo

al(position(1,j,1)) = 0.
cl(position(1,j,1)) = 2.*CX

a1(position(NX,j,l)) = 2.*AX
c1(position(NX,j,1)) = 0.

a3(position(1,j,1)) = 0.
c3(position(1,j,1)) = 2.*CZ

a3(position(NX,j,1)) = 0.
c3(position(NX,j,1)) = 2.*CZ

al (position( 1 ,j,NZ)) = 0.
cl (position(] ,j,NZ)) = 2. *CX

a1(position(NX,j,NZ)) = 2.*AX
c1 (position(NX,j,NZ)) = 0.

a3(position( 1 ,j,NZ)) = 2.*AZ
c3(position( l ,j,NZ)) = 0.

a3(position(NX,j,NZ)) = 2.*AZ
c3(position(NX,j,NZ)) = 0.

doi=2,NX-l

a2(position(i,1,1)) = 0.
c2(position(i,1,1)) = 2.*CY

a2(position(i,NY,1)) = 2.*AY
c2(position(i,NY,1)) = 0.

a3(position(i,1,l)) = 0.
c3(position(i,l,1)) = 2.*CZ

a3(position(i,NY, l )) = 0.
c3(position(i,NY,l)) = 2.*CZ

a2(position(i, 1 ,NZ)) = 0.
c2(position(i,l ,NZ)) = 2.*CY

a2(position(i,NY,NZ)) = 2.*AY
02(position(i,NY,NZ)) = 0.

101

enddo

a3(position(i, 1 ,NZ)) = 2.*AZ
c3(position(i, 1 ,NZ)) = 0.

a3(position(i,NY,NZ)) = 2.*AZ
c3(position(i,NY,NZ)) = 0.

top and bottom faces

doi=2,NX-l

enddo

do j = 2, NY-l
a3(position(i,j,1)) = 0.
c3(position(i,j,1)) = 2.*CZ

a3(position(i,j,NZ)) = 2.*AZ
03(position(i,j,NZ)) = 0.
enddo

side edges

do k = 2, NZ-l

a] (position( 1 ,1 ,k)) = 0.
c1(position(l,1,k)) = 2.*CX

a1 (position( 1 ,NY,k)) = 0.
c1(position(l,NY,k)) = 2.*CX

a1 (position(N X, l ,k)) = 2.*AX
c1(position(NX,1,k)) = 0.

a1 (position(NX,NY,k)) = 2.*AX
cl (position(NX,NY,k)) = 0.

a2(position(l,1,k)) = 0.
c2(position(1,l,k)) = 2.*CY

a2(position(l,NY,k)) = 2.*AY
c2(position(1,NY,k)) = 0.

a2(position(N X, l ,k)) = 0.
02(position(N X, 1 ,k)) = 2.*CY

a2(position(NX,NY,k)) = 2.*AY
c2(position(NX,NY,k)) = 0.

102

enddo

side faces
do j = 2, NY-l
do k = 2, NZ]
a1(position(1 ,j,k)) = 0.
c1(position(1,j,k)) = 2.*CX
a1(position(NX,j,k)) = 2.*AX
cl(position(NX,j,k)) = 0.
enddo
enddo
do i = 2, NX-l
do k = 2, NZ-l
a2(position(i,1,k)) = 0.
c2(position(i, l ,k)) = 2.*CY
a2(position(l,NY,k)) = 2.*AY
02(position(i,NY,k)) = 0.
enddo
enddo

bottom comers

a1(position(l,1 ,NZ)) = 0.
c1(position(1,1,NZ)) = 2.*CX

a1(position(NX,1,NZ)) = 2.*AX
c1(position(NX,1,NZ)) = 0.

al(position(1,NY,NZ)) = 0.
cl (position( 1 ,NY,NZ)) = 2.*CX

a1(position(NX,NY,NZ)) = 2.*AX
c1(position(NX,NY,NZ)) = 0.

a2(position(1,l ,NZ)) = 0.
c2(position(1,l ,NZ)) = 2.*CY

a2(position(NX,] ,NZ)) = 0.
02(position(NX,1,NZ)) = 2.*CY

a2(position( 1 ,NY,NZ)) = 2.*AY
02(position(1,NY,NZ)) = 0.

103

00000

a2(position(NX,NY,NZ)) = 2.*AY
02(position(NX,NY,NZ)) = 0.

a3(position(1,1,NZ)) = 2.*AZ
c3(position(l,1,NZ)) = 0.

a3(position(NX,1,NZ)) = 2.*AZ
c3(position(NX, 1 ,NZ)) = 0.

a3(position(1,NY,NZ)) = 2.*AZ
c3(position(1,NY,NZ)) = 0.

a3(position(NX,NY,NZ)) = 2.*AZ
c3(position(NX,NY,NZ)) = 0.

alter convection -- start it at t=0 (full face)

do i = 1, NX
do j = 1, NY
b3(position(i,j,l )) = BZ + Bind
b3(position(i,j,NZ)) = BZ + bBind
enddo
enddo

reorder a,b,c vectors for the x and y directions so that they are
in the proper order to be used by tridag.

call reorder(a1,b1,c1,l )
call reorder(a2,b2,c2,2)
initialize arrays and counters, switches

T = 0.
tone = 0.
ttwo = O.

COUNT = l
fluxswitch = 1.

THIS IS THE TIME LOOP.

A, B, and C vectors have been constructed for all three directions.
now we march through time implicitly, solving for the scaled
temperature field at each time step so that we may compare the
analytical solution with the experimental results taken from our data.

104

do while (COUNT.ne.0)
if(COUNT.eq.(FLASH+1)) fluxswitch=0.
STEP = 1
first time jump
call rfill(r,STEP,fluxswitch)
call flipr(r, 1)
call tridag(al ,b l ,c 1 ,r,u,N TOTAL)
call invpos(u,tone, 1)
second time jump
call rfill(r,STEP,fluxswitch)
call flipr(r,2)
call tridag(a2,b2,02,r,u,NTOTAL)
call invpos(u,ttwo,2)
third time jump-- this completes one full time step
call rfi11(r,STEP,fluxswitch)
call tridag(a3,b3,c3,r,u,NTOTAL)
call invpos(u,T,3)
output temperature at points of interest to necessary array
call output(count)
COUNT = COUNT + 1

if(COUNT.eq.(FLASH+1)) fluxswitch=0.

105

fourth time jump
call rfill(r,STEP,fluxswitch)
call flipr(r,2)
call tridag(a2,b2,c2,r,u,NTOTAL)
call invpos(u,tone,2)
fifth time jump
call rfill(r,STEP,fluxswitch)
call tridag(a3,b3,c3,r,u,NTOTAL)
call invpos(u,ttwo,3)
sixth time jump-- this completes two full time steps
call rfill(r,STEP,fluxswitch)
call flipr(r, 1)
call tridag(al,b1,cl,r,u,NTOTAL)
call invpos(u,T, 1)
outputtemperature at points of interest to necessary array
call output(count)
COUNT = COUNT + 1
if(COUNT.eq.(FLASH+1)) fluxswitch=0.
seventh time jump
call rfill(r,STEP,fluxswitch)
call tridag(a3,b3,c3,r,u,NTOTAL)

call invpos(u,tone,3)

106

c eighth time jump
call rfill(r,STEP,fluxswitch)
call flipr(r, 1)
call tridag(a1,b1,cl ,r,u,NTOTAL)
call invpos(u,ttwo, 1)
c ninth time jump» this completes three full time steps
call rfill(r,STEP,fluxswitch)
call flipr(r,2)
call tridag(a2,b2,c2,r,u,NTOTAL)
call invpos(u,T,Z)
c output temperature at points of interest to necessary array
call output(count)
' COUNT = COUNT + l

if((COUNT-l).ge.NPOINTS) COUNT = O

c this enddo loops back for the next full time step
enddo

c output temperature profiles to a file

c this section is used to get dimensional temperatures.

doi=1,NPOINTS

do j = 1, 5
Tout(i,j) = Tout(i,j)*((q*flashtime*1000.)/
+ (pcspec*zlength))
enddo
enddo
C this will put the output in a useful form for the marquardt program.

do i = 1, NPOINTS

107

00

0000000000

273

j = i+NPOINTS
k = i+2*NPOINTS

center TC
tmod(i) = Tout(i,2)
right now, set to old X TC

tmod(i) = Tout(i,3)
right now, set to old Y TC

tmod(k) = Tout(i,4)
enddo

open file for model output (scaled)
open (unit=92, file="Tmode1.txt", status = "unknown")
noise = 0.01

if (noiseswitch.eq. 1) then
call SYSTEM_CLOCK(idum)
if (idum.eq.0) pause 'idum = 0; restart'
idum = 4
write(*,*) 'idum = ‘,idum
do i = 1,NPOINTS

do j = 1,5
Tout(i,j) = Tout(i,j) + noise*gasdev(idum)
enddo
enddo
confinue
else
pause 'error with noiseswitch'
endif

do count = l, NPOINTS
time = COUNT*dt*flashtime
write(92,273) time,Tout(count, 1 ),Tout(count,2),Tout(count,3),
+ Tout(count,4),Tout(count,5)
enddo

format(1X, F7.3, 2X, 5F8.3)
close(92)

return
end

108

000

************************************
**********Sprr0gramS***************
************************************

integer function position(i,j,k)

computes the position number of a 3-space node in a single
vector array. Needed for compact matrix manipulation in tridag.

common /POSIT/ PNX,PNY,PNZ
integer i,j,k,PNX,PNY,PNZ

position = ((1-1 )*(PNY*PNZ))+((j- 1 )*PNZ)+k

return
end function

**************

subroutine reorder(a,b,c,dir)

this subroutine reorders the a, b, and c vectors for the tridag
routine to work properly.

CHANGE model size here

integer N MAX,N X,N Y,N Z
parameter (NX = 41)
parameter (NY = 41)
parameter (NZ = 41)
parameter (NMAX = 68921)

integer i,j,k,m,dir

real(8) a(NMAX),b(NMAX),c(NMAX)

real(8) tempa(NX,NY,NZ),tempb(NX,NY,NZ),tempc(NX,NY,NZ)
m = 1

call invpos(a,tempa,3)

call invpos(b,tempb,3)

call invpos(c,tempc,3)

if(dir .eq. 1) then

dok=1,NZ
doj=1,NY

109

do i = 1, NX
a(m) = tempa(i,j,k)
b(m) = tempb(i,j,k)
c(m) = tempc(i,j,k)

=m+1
enddo
enddo
enddo

elseif(dir .eq. 2) then

do k = 1, NZ
do i = 1, NX
do j = 1, NY
a(m) = tempa(i,j,k)
b(m) = tempb(i,j,k)
c(m) = tempc(i,j,k)
m=m+l
enddo
enddo
enddo

endif

return
END

**************

subroutine flipr(r,dir)

this subroutine reorders the r vector for the tridag
routine to work properly.

CHANGE model size here

integer NMAX,NX,NY,NZ
parameter (NX = 41)
parameter (NY = 41)
parameter (NZ = 41)
parameter (NMAX = 68921)

integer i,j,k,m,dir

real(8) r(NMAX)
real(8) tempr(NX,NY,NZ)

110

l
h, :5 :31: gym-"J

 

 

 

m = 1
call ianos(r,tempr,3)

if(dir .eq. 1) then

do k = 1, NZ
doj = 1, NY
do i = l, NX
r(m) = tempr(i,j,k)
m=m+1
enddo
enddo
enddo

elseif(dir .eq. 2) then

do k = 1, NZ
do i = l, NX
doj = 1, NY
r(m) = tempr(i,j,k)
m=m+1
enddo
enddo
enddo
endif
return
END
**************

subroutine invpos(in,out,dir)

this transforms the vector that is outputted by tridag into an array.
Should save on function calls.

CHANGE model size here

integer NMAX,NX,NY,NZ
parameter (NX = 41)
parameter (NY = 41)
parameter (NZ = 41)
parameter (NMAX = 68921)

111

integer i,j,k,LOCCOUNT,dir
real(8) in(NMAX),out(NX,NY,NZ)

LOCCOUNT = 0

if(dir.eq. 1) then
do k = 1, NZ
do j = 1, NY
do i = l, NX
LOCCOUNT = LOCCOUNT + 1
out(i,j,k) = in(LOCCOUNT)
enddo
enddo
enddo

elseif(dir.eq.2) then
do k = 1, NZ
do i = 1, NX
do j = 1, NY
LOCCOUNT = LOCCOUNT + 1
out(i,j,k) = in(LOCCOUNT)
enddo
enddo
enddo

elseif(dir.eq.3) then
do i = l, NX
do j = 1, NY
do k = 1, NZ
LOCCOUNT = LOCCOUNT + l
out(i,j,k) = in(LOCCOUNT)
enddo
enddo
enddo

endif

return
END

****************************

subroutine output(count)

This routine sucks the points of interest from the temperature

112

0 0 ()(5 0

vector (T) at each time step and saves it to a handy vector to
send back to the main program. This is what gets compared to the
experimental data for the Marquardt optimization.

write as an array here, output to file at end of routine
CHANGE model size here

integer NMAX,NX,NY,NZ
parameter (NX = 41)
parameter (NY = 41)
parameter (NZ = 41)
parameter (NMAX = 68921)

integer x(5),y(5),z(5),count,XTC1,YTC1,XTC2,YTC2
real(8) T(NX,NY,NZ), Tout(50000,5),xcorr,ycorr
real(8) tone(NX,NY,NZ),ttwo(NX,NY,NZ)

common IT EMPERATURE/ T,tone,ttwo,Tout
common IOTPT/ x,y,z,XTC 1 ,XTC2,YTC l ,YTC2,xcorr,ycorr

Tout(count,1) = T(x(1), y(l), z(1))
Tout(count,2) = T(x(2), y(2), 2(2))

Tout(count,3) = xcorr*T(XTC1,y(3),z(3)) +
+ (l .-xcorr)*T(XTC2,y(3),z(3))

Tout(count,4) = ycorr*T(x(4),YTCl,z(4)) +
+ (l .-ycorr)*T(x(4),YTC2,z(4))

Tout(count,5) = T(x(S), y(5), 2(5))

return
END

**************

subroutine tridag(a,b,c,r,u,n)

Thomas Algorithm for solving a tridiagonal matrix.

This is culled from Numerical Recipes in FORTRAN.
Solves for a vector u(1:n) of length n a tridiagonal linear set.
a(l :n), b(l :n), c(1:n), and r(l :n) are input vectors and are not
modified. Parameter: NMAX is the maximum expected value of n.

integer N MAX
parameter (NMAX = 100000)

113

integer n

real(8) a(n),b(n),c(n),r(n),u(n)
integer j

real(8) bet,gam(NMAX)

if(b(l).eq.0.)pause 'tridag: rewrite equations'

bet=b( 1 )

u(1 )=r(1 )/bet

do j=2,n
gam(i)=c(j-1 )/bet
bet=b0)-a0)*gam0)
if(bet.eq.0.)pause 'tridag failed (bet=0)’
u0)=(r0)-a0)*uo-1»/bet

enddo

do j=n-l,1,-l
ud)=u0)-gamo+1)*u0+1>

enddo

return

END

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

0000

00

subroutine rfill(r,STEP,fs)

This routine will fill the r-vector for use in tridag. It will

be called once during each step. Outputs the vector r, and
increments the STEP count by one. This takes the bulk of the
program time, I think, and may be worth optimizing.

NOTE: going from 3-space to vector space would save computational
time, but would also be a brain-buster.

CHANGE model size here

integer NMAX,NX,NY,NZ
parameter (NX = 41)
parameter (NY = 41)
parameter (NZ = 41)
parameter (NMAX = 68921)

integer i,j,k,STEP,position,QI,QJ

real(8) r(NMAX),Tout(50000,5)

real(8) T(NX,NY,NZ),tone(NX,NY,NZ)

real(8) ttwo(NX,NY,NZ)

real(8) Fo,Gl ,G2,Biot,q,Bind,qnd,nBind,dx,dy,dz

114

real(8) fs,QIB,QJB,bBind,anind
real(8) RA(9),RB(9),RC(9),RD(9),RE(9),RF(9),RG(9)
real(8) RH(9),RI(9),RJ(9),RK(9),RL(9),RM(9),RN(9)
real(8) RO(9),RP(9),RQ(9),RR(9)

common ITEMPERATURE/ T,tone,ttwo,Tout

common lRVECT/ RA,RB,RC,RD,RE,RP,RG,RH,RI,RJ,RK,RL,RM,RN,RO
common /RVECT2/ RP,RQ,RR,bBind,anind

common IPROPS/ Fo,G1,G2,Biot,q,qnd,Bind,nBind

common /RFU QI,QJ,QIB,QJB

common /GEOM/ dx,dy,dz

internal computations

do i = 2, NX-l
do j = 2, NY-l
do k = 2, NZ-l
r(position(i,j,k))=RA(STEP)*T(i,j,k)+RB(STEP)*
(T(i+1,j,k)+T(i-1,j,k))+RC(STEP)*(T(i,j+1,k)+
T(i,j- l,k))+RD(STEP)*(T(i,j,k+1)+T(i,j,k-l ))+
RE(STEP)*tone(i,j,k)+RF(STEP)*(tone(i+1 ,j,k)+tone
(i-1,j,k))+RG(STEP)*(tone(i,j+1,k)+tone(i,j-1,k))+
RH(STEP)*(tone(i, j,k+1 )+tone(i, j ,k- 1 ))+RI(STEP)*
ttwo(i,j,k)+RJ(STEP)*(ttwo(i+1,j,k)+ttwo(i-1,j,k))
+RK(STEP)*(ttwo(i,j+1 ,k)+ttwo(i,j- l ,k))+RL(STEP)*
(ttwo(i,j,k+1)+ttwo(i,j,k-1))
enddo

++++++++

2 face computations

r(position(i,j, 1 ))=RA(STEP)*T(i,j, l )+RB(STEP)*
(T(i+1 ,j,l)+T(i-1,j,1))+RC(STEP)*(T(i,j+l,l)+
T(i,j-l,1))+RD(STEP)*(2.*T(i,j,2))+
RE(STEP)*tone(i, j, l )+RF(STEP)*(tone(i+1,j,1)+tone
(i-l,j,1))+RG(STEP)*(tone(i,j+1 ,l )+tone(i,j-1,1 ))+
RH(STEP)*(2.*tone(i,j,2))+RI(STEP)*
ttwo(i,j,1)+RJ(STEP)*(ttwo(i+1,j,1)+ttwo(i-1,j,] ))
+RK(STEP)*(ttwo(i,j+l ,1 )+ttwo(i,j-1 ,1 ))+RL(STEP)*
(2.*ttwo(i,j,2))+RM(STEP)*(T(i,j,1))+
RN(STEP)*(tone(i,j,1))+RO(STEP)*(ttwo(i,j,1))

+++++++++

r(position(i,j,NZ)):RA(STEP)*T(i,j,NZ)+RB(STEP)*
(T(i+1,j,NZ)+T(i-1,j,NZ))+RC(STEP)*(T(i,j+l ,NZ)+
T(i,j-l,NZ))+RD(STEP)*(2.*T(i,j,NZ-1))+
RE(STEP)*tone(i,j,NZ)+RF(STEP)*(tone(i+1 ,j,NZ)+tone
(i-l,j,NZ))+RG(STEP)*(tone(i,j+1,NZ)+tone(i,j-1,NZ))+
RH(STEP)*(2.*tone(i,j,NZ-1))+RI(STEP)*

+++++

115

++++

ttwo(i,j,NZ)+RJ(STEP)*(ttwo(i+1 ,j,NZ)+ttwo(i-1 ,j,NZ))

+RK(STEP)*(ttwo(i,j+l,NZ)+ttwo(i,j-1,NZ))+RL(STEP)*

(2.*ttwo(i,j,NZ-1))+RP(STEP)*(T(i,j,NZ))+

RQ(STEP)*(tone(i,j,NZ))+RR(STEP)*(ttwo(i,j,NZ))
enddo

compute four edges

+++++++++ +++++++++

+++++++++

++++++

r(position(i,1,1))=RA(STEP)*T(i,1,1)+RB(STEP)*
(T(i+1,1,1)+T(i-1,1,1))+RC(STEP)*(2.*T(i,2,1))-1-
RD(STEP)*(2.*T(i,1,2))+
RE(STEP)*tone(i,l,1)+RF(STEP)*(tone(i+1,1 ,1 )+tone
(i-1,1,1))+RG(STEP)*(2.*tone(i,2,1))+
RH(STEP)*(2.*tone(i, 1 ,2))+RI(STEP)*
ttwo(i,1,1)+RJ(STEP)*(ttwo(i+1,1,1)+ttwo(i-1 ,1 ,1 ))
+RK(STEP)*(2.*ttwo(i,2,1))+RL(STEP)*
(2.*ttwo(i,1,2))+RM(STEP)*(T(i,1,1))+
RN(STEP)*(tone(i,],1))+RO(STEP)*(ttwo(i,1,1))

r(position(i, 1 ,NZ))=RA(STEP)*T(i, l ,NZ)+RB(STEP)*
(T(i+l,l ,NZ)+T(i-1,1,NZ))+RC(STEP)*(2.*T(i,2,NZ))
+RD(STEP)*(2.*T(i,1,NZ-l ))+

RE(STEP)*tone(i,1 ,NZ)+RF(STEP)*(tone(i+1 ,1 ,NZ)+tone
(i-1,1,NZ))+RG(STEP)*(2.*tone(i,2,NZ))+
RH(STEP)*(2.*tone(i,l,NZ-1))+RI(STEP)*
ttwo(i,],NZ)+RJ(STEP)*(ttwo(i+1,1,NZ)+ttwo(i-1 ,1 ,NZ))
+RK(STEP)*(2.*ttwo(i,2,NZ))+RL(STEP)*
(2.*ttwo(i,1,NZ—1))+RP(STEP)*(T(i,l ,NZ))+
RQ(STEP)*(tone(i,1 ,NZ))+RR(STEP)*(ttwo(i,1 ,NZ))

r(position(i,NY,l ))=RA(STEP)*T(i,NY,1)+RB(STEP)*
(T(i+l ,NY,] )+T(i-1,NY,1))+RC(STEP)*

(2.*T(i,NY-1,1 ))+RD(STEP)*(2.*T(i,NY,2))+
RE(STEP)*tone(i,NY,1)+RF(STEP)*(tone(i+l ,NY,] )+tone
(i-1,NY,]))+RG(STEP)*(2.*tone(i,NY-1,l))+
RH(STEP)*(2.*tone(i,NY,2))+RI(STEP)*
ttwo(i,NY,])+RJ(STEP)*(ttwo(i+1,NY,1)+ttwo(i-1,NY,1))
+RK(STEP)*(2.*ttwo(i,NY-1,1))+RL(STEP)*
(2.*ttwo(i,NY,2))+RM(STEP)*(T(i,N Y, 1 ))+
RN(STEP)*(tone(i,NY,1))+RO(STEP)*(ttwo(i,NY,1))

r(position(i,NY,NZ))=RA(STEP)*T(i,NY,NZ)+RB(STEP)*
(T(i+l ,NY,NZ)+T(i- l ,NY,NZ))+RC(STEP)*

(2. *T(i,NY- 1 ,NZ))+RD(STEP)*(2.*T(i,NY,NZ- 1 ))+
RE(STEP)*tone(i,NY,NZ)+RF(STEP)*(tone(i+1,NY,NZ)+tone
(i-l,NY,NZ))+RG(STEP)*(2.*tone(i,NY-1,NZ))+
RH(STEP)*(2.*tone(i,NY,NZ—1))+RI(STEP)*
ttwo(i,NY,NZ)+RJ(STEP)*(ttwo(i+1,NY,NZ)+ttwo(i-1,NY,NZ))

116

+
+
+

enddo

+RK(STEP)*(2.*ttwo(i,NY-1,NZ))+RL(STEP)*
(2.*ttwo(i,NY,NZ-1 ))+RP(STEP)*(T(1,NY,NZ))-I-
RQ(STEP)*(tone(i,NY,NZ))+RR(STEP)*(ttwo(i,NY,NZ))

two more faces

doj =2, NY-l

++++++++

++++++++

do k = 2, NZ—l
r(position( 1 ,j,k))=RA(STEP)*T( 1 ,j,k)+RB(STEP)*
(2.*T(2,j,k))+RC(STEP)*(T( 1 ,j+1 ,k)+
T(1,j-1,k))+RD(STEP)*(T(1 ,j,k+1)+T(1,j,k-1))+
RE(STEP)*tone( 1 ,j,k)+RF(STEP)*(2.*tone(2,j,k))
+RG(STEP)*(tone(l,j+1,k)+tone(1,j-1,k))+
RH(STEP)*(tone(1,j,k+l)+tone(1 ,j,k- l ))+RI(STEP)*
ttwo(l,j,k)+RJ(STEP)*(2.*ttwo(2,j,k))
+RK(STEP)*(ttwo(1,j+1,k)+ttwo(1 ,j- 1 ,k))+RL(STEP)*
(ttwo(l ,j,k+l)+ttwo(1,j,k-1 ))

r(position(NX,j,k))=RA(STEP)*T(NX,j,k)+RB(STEP)*
(2.*T(NX-l,j,k))+RC(STEP)*(T(NX,j+1,k)+
T(NX,j-1,k))+RD(STEP)*(T(NX,j,k+1)+T(NX,j,k-1))+
RE(STEP)*tone(NX,j,k)+RF(STEP)*(2.*tone
(NX-1,j,k))+RG(STEP)*(tone(NX,j+l ,k)+tone(NX,j-l,k))+
RH(STEP)*(tone(NX,j,k+1)+tone(NX,j,k-l ))+RI(STEP)*
ttwo(NX,j,k)+RJ(STEP)*(2.*ttwo(NX-1,j,k))
+RK(STEP)*(ttwo(NX,j+1 ,k)+ttwo(NX,j-1 ,k))+RL(STEP)*
(ttwo(NX,j,k+1)+ttwo(NX,j,k-l ))

enddo

four more edges

+++++++++

+ +-+ +

r(position(l,j,1))=RA(STEP)*T(1,j,1)+RB(STEP)*
(2.*T(2,j,1 ))+RC(STEP)*(T(1 ,j+1,1)+
T(1,j-1,l))+RD(STEP)*(2.*T(1,j,2))+
RE(STEP)*tone(l,j,1)+RF(STEP)*(2.*tone(2,j,1))
+RG(STEP)*(tone(l,j+l,1)+tone(1,j-1,1))+
RH(STEP)*(2.*tone(1 ,j,2))+RI(STEP)*
ttwo(l,j,1)+RJ(STEP)*(2.*ttwo(2,j,1))
+RK(STEP)*(ttwo(1 ,j+1 ,1)+ttwo(l ,j-l ,1 ))+RL(STEP)*
(2.*ttwo(1,j,2))+RM(STEP)*(T(1,j,1))+
RN(STEP)*(tone(l,j,1))+RO(STEP)*(ttwo(1,j,1))

r(position(],j,NZ))=RA(STEP)*T(1,j,NZ)+RB(STEP)*
(2.*T(2,j,NZ))+RC(STEP)*(T(1,j+1 ,NZ)+
T(1,j-1,NZ))+RD(STEP)*(2.*T(1,j,NZ-1))+
RE(STEP)*tone(1,j,NZ)+RF(STEP)*(2.*tone(2,j,NZ))
+RG(STEP)*(tone(1,j+1,NZ)+tone(1,j-1,NZ))+

117

+++++

+++++++++

+++++++++

0
:3
D...
D.
O

RH(STEP)*(2.*tone(1,j,NZ-l))+RI(STEP)*

ttwo(l ,j,NZ)+RJ(STEP)*(2.*ttwo(2,j,NZ))
+RK(STEP)*(ttwo(l,j+l,NZ)+ttwo(1,j-1,NZ))+RL(STEP)*
(2.*ttwo(1,j,NZ-1))+RP(STEP)*(T(1,j,NZ))+
RQ(STEP)*(tone(1,j,NZ))+RR(STEP)*(ttwo(1,j,NZ))

r(position(NX, j, 1 ))=RA(STEP)*T(N X,j, l )+RB(STEP)*
(2.*T(NX-l,j,l))+RC(STEP)*(T(NX,j+1,1)+
T(NX,j-1,1))+RD(STEP)*(2.*T(NX,j,2))+
RE(STEP)*tone(NX,j,1)+RF(STEP)*
(2.*tone(NX-1,j,1))+RG(STEP)*(tone(NX,j+1,l)+tone(NX,j-1 ,1 ))+
RH(STEP)*(2.*tone(NX,j,2))+RI(STEP)*
ttwo(NX,j,])+RJ(STEP)*(2.*ttwo(NX-1,j,1))
+RK(STEP)*(ttwo(NX,j+l ,1 )+ttwo(NX,j-1,1 ))+RL(STEP)*
(2.*ttwo(NX,j,2))+RM(STEP)*(T(NX,j,1))+

RN (STEP)*(tone(NX,j, 1 ))+RO(STEP)*(ttwo(N X,j, l ))

r(position(NX,j,NZ))=RA(STEP)*T(NX,j,NZ)+RB(STEP)*
(2.*T(NX-1,j,NZ))+RC(STEP)*(T(NX,j+1,NZ)+
T(NX,j-l,NZ))+RD(STEP)*(2.*T(NX,j,NZ-l))+
RE(STEP)*tone(NX,j,NZ)+RF(STEP)*
(2.*tone(NX-1,j,NZ))+RG(STEP)*(tone(NX,j+1,NZ)+
tone(NX,j-1,NZ))+RH(STEP)*(2.*tone(NX,j,NZ-1))+RI(STEP)*
ttwo(NX,j,NZ)+RJ(STEP)*(2.*ttwo(NX-l,j,NZ))
+RK(STEP)*(ttwo(NX,j+1,NZ)+ttwo(NX,j—l,NZ))+RL(STEP)*
(2.*ttwo(NX,j,NZ-1))+RP(STEP)*(T(NX,j,NZ))+
RQ(STEP)*(tone(NX,j,NZ))+RR(STEP)*(ttwo(NX,j,NZ))

last two faces

do k = 2, NZ-l
do i = 2, NX-l
r(position(i, 1 ,k))=RA(STEP)*T(i, 1 ,k)+RB(STEP)*
+ (T(i+l,1,k)+T(i-1,l,k))+RC(STEP)*(2.*T(i,2,k))
+ +RD(STEP)*(T(i,1,k+l)+T(i,1,k-1))+
+ RE(STEP)*tone(i, 1 ,k)+RF(STEP)*(tone(i+1 ,1 ,k)+-tone
+ (i- 1 , l ,k))+RG(STEP)*(2.*tone(i,2,k))+
+ RH(STEP)*(tone(i,l,k+l)+tone(i,1,k-1))+RI(STEP)*
+ ttwo(i,1,k)+RJ(STEP)*(ttwo(i+1 ,l,k)+ttwo(i-l,1,k))
+ +RK(STEP)*(2.*ttwo(i,2,k))+RL(STEP)*
+ (ttwo(i,l,k+1)+ttwo(i,1,k-1))
r(position(i,NY,k))=RA(STEP)*T(i,NY,k)+RB(STEP)*
+ (T(i+1,NY,k)+T(i-1,NY,k))+RC(STEP)*(2.*
+ T(i,NY-l,k))+RD(STEP)*(T(i,NY,k+1)+T(i,NY,k-1))+

118

++++++

RE(STEP)*tone(i,NY,k)+RF(STEP)*(tone(i+1,NY,k)+tone
(i-1,NY,k))+RG(STEP)*(2.*tone(i,NY-l ,k))+
RH(STEP)*(tone(i,NY,k+1)+tone(i,NY,k-1))+RI(STEP)*
ttwo(i,NY,k)+RJ(STEP)*(ttwo(i+1,NY,k)+ttwo(i-1,NY,k))
+RK(STEP)*(2.*ttwo(i,NY-1,k))+RL(STEP)*
(ttwo(i,NY,k+1)+ttwo(i,NY,k-1))

enddo

last four edges

++++++++ ++++++++

++++++++

++++++

r(position(l,l,k))=RA(STEP)*T(1,1,k)+RB(STEP)*
(2.*T(2,1,k))+RC(STEP)*(2.*T(1,2,k))+
RD(STEP)*(T(] ,1 ,k+1 )+T(] ,1,k-1 ))+
RE(STEP)*tone(1,1,k)+RF(STEP)*(2.*tone(2,1,k))
+RG(STEP)*(2.*tone(1,2,k))+
RH(STEP)*(tone(],1,k+1)+tone(1,1,k-l))+RI(STEP)*
ttwo(1,1,k)+RJ(STEP)*(2.*ttwo(2,1,k))
+RK(STEP)*(2.*ttwo(1,2,k))+RL(STEP)*
(ttwo(1,1,k+1)+ttwo(l ,1,k-1))

r(position( 1 ,NY,k))=RA(STEP)*T( l ,NY,k)+RB(STEP)*
(2.*T(2,NY,k))+RC(STEP)*
(2.*T(1,NY-1,k))+RD(STEP)*(T(1,NY,k+l)+T(1 ,NY,k-1))+
RE(STEP)*tone( l ,NY,k)+RF(STEP)*(2. *tone(2,NY,k))
+RG(STEP)*(2.*tone(1,NY-1,k))+

RH(STEP)*(tone(] ,NY,k+1)+tone(1 ,NY,k-l ))+RI(STEP)*
ttwo( 1 ,NY,k)+RJ(STEP)*(2. *ttwo(2,NY,k))
+RK(STEP)*(2.*ttwo(1,NY-l ,k))+RL(STEP)*
(ttwo(l,NY,k+1)+ttwo(l,NY,k-1))

r(position(NX, 1 ,k))=RA(STEP)*T(NX, 1 ,k)+RB(STEP)*
(2.*T(NX-1,l,k))+RC(STEP)*(2.*T(NX,2,k))
+RD(STEP)*(T(N X, 1 ,k-t- 1 )+T(NX, l ,k- l ))+
RE(STEP)*tone(NX, 1 ,k)+RF(STEP)*
(2.*tone(NX-1,1,k))+RG(STEP)*(2.*tone(NX,2,k))+
RH(STEP)*(tone(NX,1,k+1)+tone(NX, l ,k-l ))+RI(STEP)*
ttwo(NX, 1 ,k)+RJ(STEP)*(2.*ttwo(NX-l ,1 ,k))
+RK(STEP)*(2.*ttwo(NX,2,k))+RL(STEP)*

(ttwo(N X, l ,k+1 )+ttwo(NX, l ,k-1))

r(position(NX,NY,k))=RA(STEP)*T(NX,NY,k)+RB(STEP)*
(2.*T(NX-1,NY,k))+RC(STEP)*(2.*T(NX,NY-1,k))
+RD(STEP)*(T(NX, NY ,k+1)+T(NX, NY, k- 1))+
RE(STEP)*tone(NX, NY ,.k)+RF(STEP)*(2 *tone(NX- 1 ,NY ,k))
+RG(STEP)*(2.*tone(NX,NY-1,k))+

RH(STEP)*(tone(NX, NY, k+1 )+tone(NX,NY,k- 1 ))+RI(STEP)*
ttwo(NX,NY,k)+RJ(STEP)*(2.*ttwo(NX- l ,NY,k))

119

+
+

enddo

+RK(STEP)*(2. *ttwo(NX,NY- l ,k))+RL(STEP)*
(ttwo(NX,NY,k+1)+ttwo(NX,NY,k-1))

compute the eight corners

+++++++++ +++++++++

+++++++++

++++++

r(position(],1,1))=RA(STEP)*T(1,1,l )+RB(STEP)*
(2.*T(2,1,1))+RC(STEP)*(2.*T(1,2,1))+
RD(STEP)*(2.*T(1,1,2))+
RE(STEP)*tone(],l,1)+RF(STEP)*(2.*tone(2,1,1))
+RG(STEP)*(2.*tone(l ,2,1))+

RH(STEP)*(2.*tone( 1,1,2))+RI(STEP)*
ttwo(l,1,1)+RJ(STEP)*(2.*ttwo(2,1,1))
+RK(STEP)*(2.*ttwo(1,2,1))+RL(STEP)*
(2.*ttwo(1,1,2))+RM(STEP)*(T(1,1,1))+
RN(STEP)*(tone(],l,1))+RO(STEP)*(ttwo(1,1,1))

r(position(l,l ,NZ))=RA(STEP)*T(1,1,NZ)+RB(STEP)*
(2.*T(2,1,NZ))+RC(STEP)*(2.*T(1,2,NZ))
+RD(STEP)*(2.*T(1,1,NZ-1))+
RE(STEP)*tone(],1,NZ)+RF(STEP)*(2.*tone(2,1,NZ))+
RG(STEP)*(2.*tone(1,2,NZ))+

RH(STEP)*(2.*tone(] ,1,NZ-1))+RI(STEP)*
ttwo(l,1,NZ)+RJ(STEP)*(2.*ttwo(2,1,NZ))
+RK(STEP)*(2.*ttwo(1 ,2,NZ))+RL(STEP)*
(2.*ttwo(l,l,NZ-1))+RP(STEP)*(T(1 ,l ,NZ))+
RQ(STEP)*(tone( l , 1 ,NZ))+RR(STEP)*(ttwo(1 ,1 ,NZ))

r(position(],NY,]))=RA(STEP)*T(1,NY,1)+RB(STEP)*
(2.*T(2,NY,1))+RC(STEP)*
(2.*T(1,NY-1,1))+RD(STEP)*(2.*T(1,NY,2))+
RE(STEP)*tone( 1 ,NY, 1 )+RF(STEP)*(2. *tone(2,NY, 1 ))
+RG(STEP)*(2.*tone(1,NY-1,1))+

RH(STEP)*(2.*tone( 1 ,NY,2))+RI(STEP)*

ttwo(1,NY,1 )+RJ(STEP)*(2.*ttwo(2,NY,1))
+RK(STEP)*(2.*ttwo(1 ,NY-1,1 ))+RL(STEP)*
(2.*ttwo(1,NY,2))+RM(STEP)*(T(1,NY,1))+

RN (STEP)*(tone( 1 ,NY, 1 ))+RO(STEP)*(ttwo( 1 ,N Y, 1 ))

r(position(] ,NY,NZ))=RA(STEP)*T(1 ,NY,NZ)+RB(STEP)*
(2.*T(2,NY,NZ))+RC(STEP)*
(2.*T(l,NY-1,NZ))+RD(STEP)*(2.*T(1,NY,NZ—l))+
RE(STEP)*tone( 1 ,NY,NZ)+RF(STEP)*(2.*tone(2,NY,NZ))
+RG(STEP)*(2.*tone(1,NY-1,NZ))+
RH(STEP)*(2.*tone(1,NY,NZ-1))+RI(STEP)*

ttwo( 1 ,NY,NZ)+RJ(STEP)*(2.*ttwo(2,NY,NZ))

120

++

+++++++++ +++++++++ +++++++++

++++++++

+RK(STEP)*(2.*ttwo(1 ,NY-1 ,NZ))+RL(STEP)*
(2.*ttwo(l ,NY,NZ-l ))+RP(STEP)*(T(1,NY,NZ))+
RQ(STEP)*(t0ne(1,NY,NZ))+RR(STEP)*(ttwo( 1 ,NY,NZ))

r(position(NX,],l))=RA(STEP)*T(NX,1,1)+RB(STEP)*
(2.*T(NX-1,l,1))+RC(STEP)*(2.*T(NX,2,1))
+RD(STEP)*(2.*T(NX,1,2))+
RE(STEP)*tone(NX,1,1)+RF(STEP)*(2.*tone

(NX-1 ,1 ,1))+RG(STEP)*(2.*tone(NX,2,1))+
RH(STEP)*(2.*tone(NX, l ,2))+RI(STEP)*
ttwo(NX,l,1)+RJ(STEP)*(2.*ttwo(NX-1,l,1))
+RK(STEP)*(2.*ttwo(NX,2,1))+RL(STEP)*
(2.*ttwo(NX,l ,2))+RM(STEP)*(T(NX,1,1))+

RN (STEP)*(tone(NX, 1 ,1))+RO(STEP)*(ttwo(NX,1,1))

r(position(NX,1,NZ))=RA(STEP)*T(NX,l ,NZ)+RB(STEP)*
(2. *T(NX- l , l ,NZ))+RC(STEP)*(2.*T(NX,2,NZ))
+RD(STEP)*(2.*T(N X, 1 ,NZ- 1 ))+

RE(STEP)*tone(N X, 1 ,NZ)+RF(STEP)*(2. *tone
(NX-l,1,NZ))+RG(STEP)*(2.*tone(NX,2,NZ))+
RH(STEP)*(2.*tone(N X, 1 ,NZ- 1 ))+RI(STEP)*

ttwo(NX, 1 ,NZ)+RJ(STEP)*(2.*ttwo(NX-1,1 ,NZ))
+RK(STEP)*(2.*ttwo(NX,2,NZ))+RL(STEP)*

(2.*ttwo(N X, 1 ,NZ- 1 ))+RP(STEP)*(T(NX, l ,N Z))+
RQ(STEP)*(tone(N X, 1 ,NZ))+RR(STEP)*(ttwo(NX, 1 ,NZ))

r(position(NX,NY,1))=RA(STEP)*T(NX,NY,1)+RB(STEP)*
(2.*T(NX-1,NY,1))+RC(STEP)*
(2.*T(NX,NY-1,l))+RD(STEP)*(2.*T(NX,NY,2))+
RE(STEP)*tone(NX,NY, 1 )+RF(STEP)*(2. *tone
(NX-1,NY,]))+RG(STEP)*(2.*tone(NX,NY-1,1))+
RH(STEP)*(2.*tone(NX,NY,2))+RI(STEP)*

ttwo(NX,NY, 1 )+RJ (STEP)*(2. *ttwo(NX- 1 ,N Y, 1 ))
+RK(STEP)*(2.*ttwo(NX,NY-l ,1))+RL(STEP)*

(2. *ttwo(NX,N Y,2))+RM(STEP)*(T(NX,N Y, 1 ))+

RN (STEP)*(tone(NX,N Y, 1 ))+RO(STEP)*(ttwo(N X,NY, 1 ))

r(position(NX,NY,NZ))=RA(STEP)*T(NX,NY,NZ)+RB(STEP)*
(2.*T(NX-1,NY,NZ))+RC(STEP)*

(2.*T(NX,NY-1 ,NZ))+RD(STEP)*(2.*T(NX,NY,NZ-l ))+
RE(STEP)*tone(NX,NY,NZ)+RF(STEP)*(2.*tone
(NX-1,NY,NZ))+RG(STEP)*(2.*tone(NX,NY-1,NZ))+
RH(STEP)*(Z.*tone(NX,NY,NZ-1))+RI(STEP)*
ttwo(NX,NY,NZ)+RJ(STEP)*(2.*ttwo(NX-1,NY,NZ))
+RK(STEP)*(2.*ttwo(NX,NY-1,NZ))+RL(STEP)*
(2.*ttwo(NX,NY,NZ-1))+RP(STEP)*(T(NX,NY,NZ))+

121

0 (If? 0 (DC? 0 (>15 0 ()(7 0 ()(1 0 ()(5 0 ()(fi 0 ()(5 0 ()(5 0 (3(3 0 ()(3(1 0

+ RQ(STEP)*(tone(NX,NY,NZ))+RR(STEP)*(ttwo(NX,NY,NZ))

FLUX SECTION -- here, the illuminated area is set

do i = 1, Q1
do j = 1, Q]
r(position(i,j, 1)) = r(position(i,j, l ))+fs*qnd
enddo
enddo

***********************

DIFFERENT APERTURE WIDTHS

NOTE: this is set for a 1 1.14mm width square beam with comers
that have a ~3.17mm fillet. (5.57 xapt)

these terms add a portion of the total flux to the boundary of

the heated surface -- needed if the actual geometry doesn't line
up with the descretized geometry.

borderline edge nodes -- x-edge

do j = 1, 7
do j = 1, Q]

r(position(QI+l ,j,1)) = r(position(QI+1,j,1))+fs*QIB *qnd
enddo

borderline edge nodes -- y-edge

do i = l, 7
do i = 1, Q1

r(position(i,QJ +1 , l )) = r(position(i,QJ+1,1))+fs*QJB*qnd
enddo

borderline comer node
r(position(QI+1 ,QJ +1 , 1 )) = r(position(QI+l,QJ+1,1))+fs*QIB *QJ B *qnd
subtract out the interior nodes that are not to be heated

r(position(l 1,10,1)) = r(position(l 1,10,1))-fs*qnd
r(position( 1 1,1 1,1 )) = r(position(] 1,1 1,1 ))-fs*qnd
r(position(10,1 1.1)) = r(position( 10,1 1,1))-fs*qnd

*************************

this section for an 11.60mm width beam with comers that have a
~3.17mm fillet XAPT = 5.85

122

0 ()(5 0 ()(5 0 ()(5 0 ()(5 0 (>(5 0 ()(1 0 (3(7 0

borderline edge nodes -- x-edge

do j = l, 6
r(position(QI+1,j,l )) = r(position(QI+1,j,1))+fs*QIB *qnd
enddo

borderline edge nodes -- y-edge

do i = 1, 6
r(position(i,QJ +1 , 1 )) = r(position(i,QJ+],1))+fs*QJB*qnd
enddo

subtract out the interior nodes that are not to be heated

r(position(12,10,1)) = r(position(12,10,l ))-fs*qnd
r(position(12,1 l ,1)) = r(position(12,1 1,1))-fs*qnd
r(position(12,12,l )) = r(position(12,12,1))-fs*qnd
r(position(10,12,1)) = r(position(10,12,1))-fs*qnd
r(position(] 1,12,1)) = r(position(l 1,12,1))-fs*qnd

r(position(] 1,1 1,1)) = r(position(] 1,1 1,1))-0.5*fs*qnd
r(position(12,9,l)) = r(position(12,9,1))-0.35*fs*qnd
r(position(9,] 2,1)) = r(position(9,]2,1))-0.35*fs*qnd

********************

this section for an 10.62mm width beam with corners that have a
~3.17mm fillet (5.31 xapt)
borderline edge nodes -- x-edge

do j = 1, 4
r(position(QI+1,j,l )) = r(position(QI+1,j,l))+fs*QIB *qnd
enddo

borderline edge nodes -- y-edge

do i = 1, 4
r(position(i,QJ +1 , l )) = r(position(i,QJ+1,1))+fs*QJB*qnd
enddo

subtract out the interior nodes that are not to be heated

r(position(] 1,10,1)) = r(position(l l ,10,1))-fs*qnd
r(position(] 1,1 1, l )) = r(position(] 1,1 1 ,1))-fs*qnd
r(position(10,l 1,1)) = r(position(10,1 1,1))-fs*qnd
r(position(9,] 1,1)) = r(position(9,1 1,1 ))-fs*qnd
r(position(] 1 ,9,1)) = r(position(] 1,9,1))-fs*qnd

123

0 (DC? 0 (DC? 0 ()(5 0 (DC? 0 (it? 0 (3(3 0 ()(5 0 (3!? 0 ()(7 0 ()(7 0 ()(1 0 ()(7 0 (if? 0 ()(5 0 ()(3 0

r(position(10,10, 1)) = r(position(10,10,1))-fs*qnd

r(position(8,] 1,1)) = r(position(8,1 1,1))—0.5 *fs*qnd
r(position(] 1,8,1 )) = r(position(] l ,8,1))-0.5*fs*qnd

****************************

NOTE: this is set for an 11.30mm width square beam with comers
that have a ~3.17mm fillet. XAPT = 5.65

these terms add a portion of the total flux to the boundary of

the heated surface -- needed if the actual geometry doesn't line

up with the descretized geometry.

borderline edge nodes -- x-edge

do j = l, 6
r(position(QI+l,j,1)) = r(position(QI+1,j,1))+fs*QIB *qnd
enddo

borderline edge nodes -- y-edge

do i = 1, 6
r(position(i,QJ +1 , 1 )) = r(position(i,QJ +1 , l ))+fs*QJ B*qnd
enddo

add nodes on edge that are partially heated
r(position(12,9,l)) = r(position(12,9,1))+0.1*fs*QIB *qnd
r(position(9,]2,1)) = r(position(9,]2,1))+0.1*fs*QJB*qnd
r(position(12,8,l)) = r(position(12,8,1))+0.55 *fs*QIB*qnd
r(position(8,]2,1)) = r(position(8, 12,1 ))+0.55*fs*QJB*qnd
r(position(12,7,1)) = r(position(12,7,l ))+0.95*fs*QIB*qnd
r(position(7,12, 1 )) = r(position(7,12,1))+0.95*fs*QJB *qnd

subtract out the interior nodes that are not to be heated

r(position(] 1,1 1,1)) = r(position(l 1,1 1, l ))-fs*qnd
r(position(] 1,10,1)) = r(position(] 1,10, l ))-0.8*fs*qnd
r(position(10,1 1, l )) = r(position(10,1 1,1))-0.8*fs*qnd

********************

OLD ALTER APERTURE SIZE CODE
this is the section where nodes are suntracted off (if not merely
commented out above).

i=12

do i = 1 1,12
r(position(i, 1 3,1)) = r(position(i,13,1))-fs*QIB*qnd

124

0 (DC? 0 ()(5 0 ()(3 0 ()(5 0 ()(3 0 (>(5 0 ()(5 0 (ht? 0 (3!? 0 (3(5 0

0

r(position(13,i,1)) = r(position(13,i,1 ))-fs*QJB*qnd
enddo

r(position(12,12,1)) = r(position(12,12,1))-fs*qnd

EDWIN-LIKE APERTURE CODE
fix convection to make it like Edwin's model

do i = 1, 6
do j = 1, 6
r(position(i,j,1)) = r(position(i,j,1))
+ -fs*RM(STEP)*(T(i,j,1))—fs*RN(STEP)*(tone(i,j,1))
-fs*RO(STEP)*(ttwo(i,j, l ))
enddo
enddo
do i = 1, 6

r(position(i,7,1)) = r(position(i,7, l ))-fs*0.65 *RM(STEP)*
+ (T(i,7,1))-fs*0.65*RN(STEP)*(tone(i,7,1))-fs*0.65*RO(STEP)*
(ttwo(i,7, l ))
enddo

do j = l, 6
r(position(7,j,1)) = r(position(7,j, 1 ))-fs*0.65 *RM(STEP)*
+ (T(7,j,1))-fs*0.65 *RN(STEP)*(tone(7,j,1))-fs*0.65*RO(STEP)*
(ttwo(7,j,1))
enddo

r(position(7,7,1)) = r(position(7,7,1))-fs*0.65 *0.65 *RM(STEP)*
+(T(7,7,1))-fs*0.65 *0.65 *RN(STEP)*(tone(7,7, l ))-fs*0.65*0.65*

+RO(STEP)*(ttwo(7,7,1))

increment the time step
STEP = STEP + 1

return
END

***********************
subroutine modelflux(epsil4,t4pos,t4neg)
this routine reads the model output and adjusts it based on the

heat flux adjustments used in the sensitivity coefficient portion
of the program.

125

L _:__zs_ ~ luau-l

 

 

0 ()(5 O

real(8) tmod(50000),atry(7),epsil4,dumpos,dumneg,flashtime
real(8) t4pos(50000),t4neg(50000),zlength
integer nread,i

common IFLASHIT/ flashtime
common IMAINPROG/ tmod,atry
common IFOO/ nread

common /ZL/ zlength

dumpos = atry(4)+epsil4
dumneg = atry(4)-epsil4

re-nondimensionalize the temperature history

do i = 1, (3*nread)
tmod(i) = tmod(i)/((atry(4)*flashtime*1000)/zlength)

t4pos(i) = tmod(i)*((dumpos*flashtime*1000)/zlength)
t4neg(i) = tmod(i)*((dumneg*flashtime*1000)/zlength)
t4pos(i) = (tmod(i)*dumpos)/atry(4)
t4neg(i) = (tmod(i)*dumneg)/atry(4)

enddo

return
end

**************************

FUNCTION gasdev(idum)
integer idum

*uses ran2*
returns a normally distributed deviate with zero mean and unit

variance, using ran2(idum) as the source of uniform deviates
culled from Numerical Recipes in FORTRAN

integer iset

real(8) fac,gset,rsq,v1,v2

save iset,gset

data iset/O/

if(iset.eq.0) then
v 1 =2. *ran2(idum)-1.
v2=2. *ran2(idum)—1.
rsq=v1**2+v2**2
if(rsq.ge. 1 ..or.rsq.eq.0.) goto l
fac=sqrt(-2.*log(rsq)/rsq)
gset=v1*fac

126

11

gasdev=v2*fac
iset=1
else
gasdev=gset
iset=0
endif

return
end

316*************************

FUNCTION RAN2(IDUM)

from Numerical Recipes in FORTRAN
returns a uniform random deviate between 0. and 1.

PARAMETER (M=714025,IA=1366,IC=150889,RM=1 .4005] 12E-6)
DIMENSION IR(97)

DATA IFF /O/

IF(IDUM.LT.0.0R.IFF.EQ.O)THEN
IFF=1
IDUM=MOD(IC-IDUM,M)

DO 1 1 J =1 ,97
IDUM=MOD(IA*IDUM+IC,M)
IR(J)=IDUM
CONTINUE

IDUM=MOD(IA*IDUM+IC,M)

IY=IDUM

ENDIF

J=1+(97*IY)/M

IF(J.GT.97.0R.J.LT.1)PAUSE

IY=IR(J)

RAN2=IY*RM

IDUM=MOD(IA*IDUM+IC,M)

IR(J)=IDUM

RETURN

END

****************************

end of subprograms, file

127

APPENDIX C

RANDOM IN -PLANE AN ISOTROPY

FINITE DIFFERENCE PROGRAM

128

The following pages contain FORTRAN code for an explicit finite difference program
that models the diffusion equation and assumes monoclinic heat conduction in a slab—like

specimen. The thermal diffusivity tensor 0 that allows a monoclinic response may be
expressed as

0111 0112 0
a: 0121 0122 0 ((2.1)
O O 0.33

where 0.12 = 0121. The parameters to be estimated are the four nonzero components of a.
The simulated specimen boundary conditions consist of: insulated edges, convection on
the top and bottom surfaces (of different magnitude), and a circular (top surface) heat
pulse with a 6 mm radius of constant magnitude and finite duration. The model verifies .
that the Courant stability criterion is not violated before each run. The simulated
geometry is easily adjustable, with the exception of the geometry of the illuminated area,
where the weighting parameters must be adjusted any time the flash diameter is changed.
Most of the experimental and geometric parameters are adjusted via the setup subroutine
that precedes the actual model program. The output consists of an array containing the

temperature history of the simulated specimen at four points on the bottom surface.

129

 

 

000000

Sean Davis -- programmer

this program computes the temperature history of a single layer
specimen that is heated over a portion of the top face. The model
handles the fully populated thermal diffusivity tensor, which
includes mixed derivatives. The top surface experiences convection,

while the other surfaces are either insulated or experience symmetric
boundary conditions.

subroutine setup(IMAX,stopiter)

CHANGE model size here
integer NX,NY,NZ

parameter (NX = 81)
parameter (NY = 81)
parameter (NZ = 41)

declare variables

integer IMAX,sensors,x(5),y(5),z(5),XTC,YTC,nreada,nread,i, j
integer dataswitch,noiseswitch,MI,MJ,nn

real(8) axx,ayy,azz,axy,pcspec,qest,topconv,botconv,outtime
real(8) flashtime,fliz,stopiter,zlength,xlength,ylength,xapt,lxtc
real(8) lx,ly,dz,dx,dy,xcorr,ycorr,tdiff(7),timecalc

real(8) tdiffo(7),Tmax,noise,qn(NX,NY)

character*12 fname
common blocks

common /GUESS 1/ tdiff,tdiffo

common /OTPT/ x,y,z,XTC,YTC,xcorr,ycorr,noise,noiseswitch
common IAPT/ qn

common IREADITT/ nreada,dataswitch,fliz,nn
common /GEOM/ dx,dy,dz,outtime,pcspec
common IT [MEI timecalc

common /FOO/ nread,lxtc,xapt

common /FLASHIT/ flashtime

common IFISH/ sensors,Tmax

common /FILE/ fname

130

0000000000

000000000

00000

0

set properties -- for inverse problem, these are guesses

m1841 parameters

azz = 1.1D-07

principal directions (0 degrees)
axx = 2.62D-07

ayy = 0.58D-07

axy = 0.0D-07

max shear (45 degrees)

axx = 1.6D-07

ayy = 1.6D-07

axy = 1.02D-07

m 1842 parameters

azz = 1.1D-07
axx = 1.8D-07
ayy = 1.4D-07
axy = 0.2D-07
axx = 1.883D-07
ayy = 1.317D-07
axy = 0.0D-07

1843 params

azz = 1.1D-07
axx = 1.600D-07
ayy = 1.600D-07
axy = 0.564D-07

inverse parameters

azz = 1.00D-07

axx = 1.00D-07

ayy = 1.00D-07

axy = 0.50D-07

use axy below if setting shear to zero
axy = 0.00D-07

these properties are assumed known
pcspec = 1.54D+06
qest = 0.025

topconv = 1.3D-05
botconv = 5./pcspec

131

set experiment parameters
fname = "m1843e.txt"

timecalc = 100.
flashtime = 1.0
fhz = 30.

for data input, read every "nn"th point (nn=l means'read every datafile point)
nn = 1

outtime = 1./fhz
nreada = 3010

IMAX = 50
stopiter = 1.0D-10

sensors = 4
Tmax = 0.0

dataswitch: l = inverse problem, 2 = forward problem
dataswitch = 1

forward problem noise parameters -- switch and level
noiseswitch: l = noise, 0 = clean data
(a noise level of 0.01 = IR Camera)

noiseswitch = 0
if (dataswitch.eq.1) noiseswitch = 0
noise = 0.01

set geometry and output information

zlength = 3.25
xlength = 40.0
xapt = 6.00

hm=90

zlength = zlength/1000.
xlength = xlength/1000.
ylength = xlength

1x = lxtc/ 1000.
1y = 1x

132

dz = zlength/dfloat(NZ- 1)
dx = xlength/dfloat(NX- l)
dy = ylength/dfloat(NY- l)

compute center point of model

M1 = (NX+1)/2
M1 = (NY+1)/2

for output nodes, "0" nodes are overwritten to account for the
inter-nodal sensor locations

x(1)= MI
x(2) = 0
x(3) = M1
x(4) = 0
x(5) = M1
x(6) = 0
x(7) = 1

y(1)= MI
W) = M1
y(3) = 0
y(4) = 0
y(5) = M1
y(6) = 1
y(7) = 0

2(1) = NZ
2(2) = NZ
2(3) = NZ
z(4) = NZ
2(5) = 1
2(6) = 1
2(7) = 1

XTC = M1 + int(lx/dx)
YTC = M] + int(ly/dy)

xcorr = 1. - ((lx - float(XTC-MI)*dx)/dx)
ycorr = l. - ((ly - float(YTC-MJ)*dy)/dy)

set flux absorption constants -- this models a full circle
with a radius of 6.0mm

133

qn=0.

do j = MJ-3, MJ+3
do i = MI-l l, MI+11
qn(i,j) = 1.
enddo
enddo

do i = MI-lO, MI+10
qn(i,MJ+4) = 1.
qn(i,MJ+5) = 1.
qn(i,MJ-4) = 1.
qn(i,MJ-S) = 1.
enddo

do i = MI-9, MI+9
qn(i,MJ+6) = 1.
qn(i,MJ+7) = l.
qn(i,MJ-6) = 1.
qn(i,MJ-7) = l.

 

enddo

do i = MI-8, MI+8
qn(i,MJ+8) = 1.
qn(i,MJ-8) = 1.
enddo

do i = MI-7, MI+7
qn(i,MJ+9) = l.
qn(i,MJ-9) = 1.
enddo

do i = MI-5, MI+5
qn(i,MJ+10) = 1.
qn(i,MJ-10)= 1.
enddo

do i = MI-3, MI+3
qn(i,MJ+l 1) = 1.
qn(i,MJ-l l) = 1.
enddo

qn(MI,MJ+12) = 0.5

qn(MI+12,MJ) = 0.5
qn(MI+1 ,MJ+12) = 0.45

134

qn(MI+12,MJ+1) = 0.45
qn(MI+12,MJ+2) = 0.35
qn(MI+2,MJ+12) = 0.35
qn(MI+11,MJ+4) = 0.8
qn(MI+4,MJ+l l) = 0.8
qn(MI+6,MJ+10) = 0.85
qn(MI+10,MJ+6) = 0.85
qn(MI+10,MJ+7) = 0.2
qn(MI+7,MJ +10) = 0.2
qn(MI+9,MJ+8) = 0.5
qn(MI+8,MJ+9) = 0.5
qn(MI+12,MJ+3) = 0.1
qn(MI+3,MJ+12) = 0.1
qn(MI+11,MJ+5) = 0.4
qn(MI+5,MJ+1 1) = 0.4

qn(MI,MJ-12) = 0.5
qn(MI-12,MJ) = 0.5
qn(MI- l ,MJ-12) = 0.45
qn(MI-12,MJ-l) = 0.45
qn(MI-12,MJ-2) = 0.35
qn(MI-2,MJ-12) = 0.35
qn(MI-l 1,MJ-4) = 0.8
qn(MI-4,MJ-1 l) = 0.8
qn(MI—6,MJ -10) = 0.85
qn(MI-10,MJ-6) = 0.85
qn(MI-10,MJ-7) = 0.2
qn(MI-7,MJ-10) = 0.2
qn(MI-9,MJ-8) = 0.5
qn(MI-8,MJ-9) = 0.5
qn(MI-12,MJ-3) = 0.1
qn(MI-3,MJ-l2) = 0.1
qn(MI-l l,MJ-5) = 0.4
qn(MI-5,MJ-l 1) = 0.4

qn(MI+1,MJ-12) = 0.45
qn(MI+12,MJ-l) = 0.45
qn(MI+12,MJ-2) = 0.35
qn(MI+2,MJ-12) = 0.35
qn(MI+1 l,MJ-4) = 0.8
qn(MI+4,M.1-l 1) = 0.8
qn(MI+6,MJ-10) = 0.85
qn(MI+10,MJ-6) = 0.85
qn(MI+10,MJ-7) = 0.2
qn(MI+7,MJ-10) = 0.2
qn(MI+9,MJ-8) = 0.5

135

 

 

Onoooono

0

qn(MI+8,MJ-9) = 0.5

qn(MI+12,MJ-3) = 0.1
qn(MI+3,MJ-12) = 0.1
qn(MI+1 l,MJ-5) = 0.4
qn(MI+5,MJ-l 1) = 0.4

qn(MI—1,MJ+12) = 0.45
qn(MI-12,MJ+1) = 0.45
qn(MI-12,MJ+2) = 0.35
qn(MI-2,MJ+12) = 0.35
qn(MI-11,MJ+4) = 0.8
qn(MI-4,MJ+1 l) = 0.8
qn(MI-6,MJ+10) = 0.85
qn(MI-10,MJ+6) = 0.85
qn(MI-10,MJ+7) = 0.2
qn(MI-7,MJ+10) = 0.2
qn(MI-9,MJ+8) = 0.5
qn(MI-8,MJ+9) = 0.5
qn(MI-12,MJ+3) = 0.1
qn(MI-3,MJ+12) = 0.1
qn(MI-11,MJ+5) = 0.4
qn(MI-5,MJ+11) = 0.4

use this to get l-D flash across entire surface

doi= 1,NX
doj=1,NY
qn(i,j)=1
enddo
enddo

set Marquardt estimation parameters

tdiff( 1 )=azz
tdiff(2)=axx
tdiff(3)=ayy
tdiff(4)=axy
tdiff(5)=qest
tdiff(6)=topconv
tdiff(7)=botconv

return
end

136

 

 

0

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

subroutine model
c CHANGE model size here
integer NX,NY,NZ

parameter (NX = 81)
parameter (NY = 81)
parameter (NZ = 41)

c declare variables

c real(8) tmod(25000)
real(8) dx,dy,dz,check,qnode,qout,hout,noise,tmod(80000),xapt

real(8) T(NX,NY,NZ),Tnew(NX,NY,NZ),Tout(40000,5),timeout(40000)
real(8) axx,ayy,azz,axy,pcspec,kzz,delf,bkhout,lxtc,qn(NX,NY)

real(8) fxx,fyy,fzz,fxy,hnode,time,outtime,dtf,del,atry(7)

real(8) timecalc,flashtime,dt,fliz,xcorr,ycorr,backhnode

integer FLASH,FINAL,FINAL2,XTC,YTC,COUNT,nreada,dataswitch,nn
integer i,j,k,l,m,x(5),y(5),z(5),noiseswitch,idum,nread,nreadm

c common blocks

common /MAINPROG/ tmod,atry
common /OTPT/ x,y,z,XTC,YTC,xcorr,ycorr,noise,noiseswitch

common /APT/ qn
common /GEOM/ dx,dy,dz,outtime,pcspec
common IT IME/ timecalc

common IFLASHIT/ flashtime
common IREADI'I'T/ nreada,dataswitch,fhz,nn

common IFOO/ nread,lxtc,xapt

set the physical parameters for the problem
all lengths are in millimeters, all times are in seconds.

axx = atry(2)
ayy = atry(3)
azz = atry( l)
axy = atry(4)
kzz = azz*pcspec

qout = atry(5)

137

 

_'O-.'h fisl‘p'
s

.3. ..‘

 

 

hout = atry(6)
bkhout = atry(7)

compute/check the stability criterion for the problem
this criteria is the same as that found for the simpler orthogonal

problem.

check = 1.0
dt = 2./fl1z
m = 0

do while (check. ge.0.95)
dt = dt/2.

fxx = axx*dt/(dx*dx)
fyy = ayy*dt/(dy*dy)
fzz = azz*dt/(dz*dz)
fxy = axy*dt/(2.*dx*dy)

check = 2.*fxx + 2.*fyy + 2.*fzz
m=m+1
enddo

de1= dt/lOO.

this section sets the flash timestep

dtf = 0.02/fhz
check = 1.0
m = 0

do while (check.ge.0.95)
dtf = dtf/2.

fxx = axx*dtf/(dx*dx)
fyy = ayy*dtf/(dy*dy)
fzz = azz*dtf/(dz*dz)
fxy = axy*dtf/(2.*dx*dy)

check = 2.*fxx + 2.*fyy + 2.*fzz
m=m+1
enddo

delf = dtf/ 1 OO.

compute the value of I needed to reach the end of the flash

138

and the end of the problem

FLASH = nint(flashtime/dtf)

FINAL = nint((timecalc - (2.*flashtime))/dt)
FINAL2 = nint(timecalc/dt)

compute other dt-dependent parameters
qnode = (2.*qout*dtf)/dz

hnode = (2.*hout*dtf)/dz

backhnode = (2.*bkhout*dtf)/dz

COUNT = l

T = 0.
Tnew = 0.

TIME LOOP SHALL BEGIN HERE -- CALCULATION BEGINS
loop over experiment time

do 1 = 1, 2*FLASH

if we are in the flashtime, add the heat flux to the top nodes

if(l.le.FLASH) then
do i = 1, NX
do j = 1, NY
T(i,j,l) = T(i,j,1)+qn(i,j)*qnode
enddo
enddo
endif

here we compute the temperature field in the body

four top comers

Tnew(l,l,1) = fxx*(2.*T(2,l,1)-2.*T(l,1,l))+fyy*(2.*T(l,2,1)-
+ 2.*T(1,1,1))+fzz*(2.*T(],1,2)-2.*T(l,1,1))+T(1,l,1)-hnode*T(l,1,1)

Tnew(NX,l,1) = fxx*(2.*T(NX-l,l,1)-2.*T(NX,1,1))+fyy*(2.*
+ T(NX,2,])-2.*T(NX,1,1))+fzz*(2.*T(NX,l,2)-2.*T(NX,1,1))+T(NX,1,1)

139

 

 

+ -hnode*T(NX,l , 1)

Tnew(l,NY, l) = fxx*(2.*T(2,NY,l)-2.*T(1,NY,1))+fyy*(2.*
+ T(l,NY-l,l)-2.*T(1 ,NY,1))+fzz*(2.*T(1,NY,2)-2.*T(1,NY,1 ))+
+ T(1,NY,1 )-hnode*T(1,NY, 1)

Tnew(NX,NY, 1) = fxx*(2.*T(NX-1,NY,])-2.*T(NX,NY,1))+fyy*(2.*
+ T(NX,NY-1,1)-2.*T(NX,NY,l))+fzz*(2.*T(NX,NY,2)-2.*T(NX,NY,1))+
+ T(NX,NY, 1 )-hnode*T(NX,NY, 1)

two top edges ‘2

do i = 2, NX-l

Tnew(i,1,l) = fxx*(T(i+1,1,1)-2.*T(i,1,1)+T(i-1,1,l))+
+ fyy*(2.*T(i,2,1)-2.*T(i,l,l))+fzz*(2.*T(i,l,2)-2.*T(i,l,1))+ 1:
+ T(i, l , l )-hnode*T(i, l , 1) i1

 

Tnew(i,NY,1) = fxx*(T(i+1,NY,])-2.*T(i,NY,1)+T(i-l ,NY, 1 ))+
+ fyy*(2.*T(i,NY-l ,1)-2.*T(i,NY,1))+fzz*(2.*T(i,NY,2)-
+ 2.*T(i,NY,1))+T(i,NY,1)-hnode*T(i,NY,1)
enddo

two top edges

doj = 2, NY-l
Tnew(l,j,1) = fxx*(2.*T(2,j,1)-2.*T(1,j,1))+
+ fyy*(T(l,j+1,l)-2.*T(1,j,1)+T(l,j-1,1))+fzz*
+ (2.*T(l,j,2)-2.*T(1,j,l))+T(1,j,l)-hnode*T(l,j,l)

Tnew(NX,j,1) = fxx*(2.*T(NX-l,j,1)-2.*T(NX,j,1))+
+ fyy*(T(NX,j+1,1)-2.*T(NX,j,1)+T(NX,j-l,l))+fzz*
+ (2.*T(NX,j,2)-2.*T(NX,j,1))+T(NX,j,1)-hnode*T(NX,j,l)

enddo
top face
do i = 2, NX-l

do j = 2, NY-l

Tnew(i,j,1) = fxx*(T(i+l,j,1)-2.*T(i,j,1)+T(i-1,j,1))+

+ fyy*(T(i,j+l,l)-2.*T(i,j,l)+T(i,j-l,1))+fzz*(2.*T(i,j,2)-
+ 2.*T(i,j,1))+fxy*(T(i+l,j+1,1)—T(i+1,j-1,1)-T(i-1,j+1,1)+
+ T(i-l ,j-1,l))+T(i,j,1)-hnode*T(i,j,l)

enddo
enddo

two side faces

140

 

do i = 2, NX-l

do k = 2, NZ-l
Tnew(i, 1 ,k) = fxx*(T(i+1,l,k)-2.*T(i,1,k)+T(i-1 ,1,k))+
+ fyy*(2.*T(i,2,k)-2.*T(i,1,k))+fzz*(T(i,1,k+l)-2.*T(i,1,k)+
+ T(i,] ,k-1))+T(i, 1 ,k)
Tnew(i,NY,k) = fxx*(T(i+1,NY,k)-2.*T(i,NY,k)+T(i-1,NY,k))+
+ fyy*(2.*T(i,NY-1,k)-2.*T(i,NY,k))+fzz*(T(i,NY,k+1)-2.*
+ T(i,NY,k)+T(i,NY,k-1))+T(i,NY,k)
enddo l
enddo

1
two side faces '

 

do j = 2, NY-l
do k = 2, NZ-l
Tnew( l ,j,k) = fxx*(2.*T(2,j,k)-2.*T( 1 ,j,k))+fyy*
+ (T(1,j+1,k)-2.*T(l,j,k)+T(],j-l,k))+fzz*(T(1,j,k+1)-2.*
+ T(l ,j,k)+T(],j,k-1))+T(l,j,k)
Tnew(NX,j,k) = fxx*(2.*T(NX-1,j,k)-2.*T(NX,j,k))+fyy*
+ (T(NX,j+1,k)-2.*T(NX,j,k)+T(NX,j-l,k))+fzz*(T(NX,j,k+l )-
+ 2. *T(NX,j,k)+T(NX,j,k- 1 ))+T(NX,j,k)
enddo
enddo

four side edges

do k = 2, NZ]

Tnew(] , 1,10 = fxx*(2.*T(2,l,k)-2.*T(1,1,k))+fyy*
+ (2.*T(1,2,k)-2.*T(1,1,k))+fzz*(T(],l,k+1)-2.*
+ T(1,1,k)+T(1,l,k-1))+T(1,1,k)

Tnew(NX, l ,k) = fxx*(2.*T(NX-1,1,k)-2.*T(NX,l,k))+fyy*
+ (2. *T(NX,2,k)-2.*T(N X, l ,k))+fzz*(T(N X, l ,k+1 )-
+ 2.*T(NX,1,k)+T(NX,l,k-1))+T(NX,1,k)

Tnew( 1 ,NY,k) = fxx*(2. *T(2,NY,k)-2.*T( 1 ,NY,k))-l-
+ fyy*(2.*T(1,NY-1,k)-2.*T(1,NY,k))+fzz*(T(1,NY,k+1)-2.*
+ T(l,NY,k)+T(l,NY,k-1))+T(1,NY,k)

Tnew(NX,NY,k) = fxx*(2.*T(NX-1,NY,k)-2.*T(NX,NY,k))+
+ fyy*(2.*T(NX,NY-1,k)-2.*T(NX,NY,k))+fzz*(T(NX,NY,k+1)-2.*
+ T(NX,NY,k)+T(NX,NY,k- 1 ))+T(NX,NY,k)
enddo

141

model interior

do i = 2, NX-l
do j = 2, NY-l
do k = 2, NZ-l
Tnew(i,j,k) = fxx*(T(i+1,j,k)-2.*T(i,j,k)+T(i-1,j,k))+
+ fyy*(T(i,j+l ,k)-2.*T(i,j,k)+T(i,j-1 ,k))+fzz*
+ (T(i,j,k+1)-2.*T(i,j,k)+T(i,j,k-l ))+fxy*(T(i+l ,j+1 ,k)-
+ T(i+1,j-l,k)-T(i-l,j+1,k)+T(i-1,j-1,k))+T(i,j,k)
enddo
enddo
enddo

bottom face

do i = 2, NX-l
do j = 2, NY-l
Tnew(i,j,NZ) = fxx*(T(i+1,j,NZ)-2.*T(i,j,NZ)+T(i-1,j,NZ))+

+ fyy*(T(i,j+l ,NZ)-2.*T(i,j,NZ)+T(i,j-l ,NZ))+fzz*(2.*
+ T(i,j,NZ-l)-2.*T(i,j,NZ))+fxy*(T(i+1,j+1,NZ)-
+ T(i+1,j-l,NZ)-T(i-l,j+1,NZ)+T(i-1,j-l ,NZ))+T(i,j,NZ)-
+ backhnode*T(i,j,NZ)
enddo
enddo

two bottom edges

do i = 2, NX-l

Tnew(i,],NZ) = fxx*(T(i+l,l,NZ)-2.*T(i,],NZ)+T(i-l,1,NZ))+
+ fyy*(2.*T(i,2,NZ)-2.*T(i,l ,NZ))+fzz*(2.*T(i,l,NZ-1)-
+ 2.*T(i, 1 ,NZ))+T(i, 1 ,NZ)-backhnode*T(i, 1 ,NZ)

Tnew(i,NY,N Z) = fxx*(T(i+1,NY,NZ)-2.*T(i,NY,NZ)+T(i-1,NY,NZ))+
+ fyy*(2.*T(i,NY-l ,NZ)-2.*T(i,NY,NZ))+fzz*(2.*T(i,NY,NZ-1)-
+ 2.*T(i,NY,NZ))+T(i,NY,NZ)-backhnode*T(i,NY,NZ)
enddo

two bottom edges
do j = 2, NY-l
Tnew( l ,j,NZ) = fxx*(2.*T(2,j,NZ)-2.*T( 1 ,j,NZ))+
+ fyy*(T(1,j+1,NZ)-2.*T(l,j,NZ)+T(1,j-1,NZ))+fzz*
+ (2.*T(l,j,NZ-l)-2.*T(1,j,NZ))+T(1,j,NZ)-backhnode*T(1,j,NZ)

Tnew(NX,j,NZ) = fxx*(2.*T(NX- l ,j,NZ)-2.*T(NX,j,NZ))+

142

+
+

enddo

fyy*(T(NX,j+1 ,NZ)-2.*T(NX,j,NZ)+T(NX,j- l ,NZ))+fzz*(2.*
T(NX,j,NZ- 1 )-2. *T(NX,j,NZ))+T(N X,j,N Z)-backhnode*T(NX, j,N Z)

four bottom comers

+
+
+

Tnew( l , l ,N Z) = fxx*(2.*T(2,1,NZ)-2.*T(l ,1,NZ))+fyy*(2.*T(l,2,NZ)-
2.*T(l,1,NZ))+fzz*(2.*T(1,1,NZ-1)-2.*T(1,1,NZ))+T(1,1,NZ)-
backhnode*T(1,1 ,NZ)

Tnew(NX, 1 ,NZ) = fxx*(2.*T(NX-l,1,NZ)-2.*T(NX,1,NZ))+fyy*(2.*
T(NX,2,NZ)-2.*T(NX,l,NZ))+fzz*(2.*T(NX,1,NZ-l)-2.*T(NX,1,NZ))+
T(NX, 1 ,N Z)-backhnode*T(N X, 1 ,NZ)

Tnew( 1 ,NY,NZ) = fxx*(2.*T(2,NY,NZ)-2.*T( l ,NY,NZ))+fyy*(2. *
T( l ,NY-1 ,NZ)-2.*T(1,NY,NZ))+fzz*(2.*T(1 ,NY,NZ-l )-
2.*T(1,NY,NZ))+T(1,NY,NZ)-backhnode*T(l,NY,NZ)

Tnew(NX,NY,NZ) = fxx*(2.*T(NX-1,NY,NZ)-

2. *T(NX,NY,NZ))+fyy*(2. *T(NX,NY- 1 ,NZ)-
2.*T(NX,NY,NZ))+fzz*(2.*T(NX,NY,NZ-1)-
2.*T(NX,NY,NZ))+T(NX,NY,NZ)-backhnode*T(N X,NY,N Z)

send Tnew to T once calculations are complete

T = Tnew

if the calculated time equals the desired output time, write the
temperature at the points of interest to a separate array

time = dtf*l
outtime = (1./fhz)*COUNT

if((time+delf).ge.outtime) then

timeout(COUNT) = time
Tout(COUNT,1) = T(x(1), y(l), z(1))

Tout(COUNT,2) = xcorr*T(XTC,y(2),z(2)) +
(1.-xcorr)*T((XTC+1),y(2),z(2))

Tout(COUNT,3) = ycorr*T(x(3),YTC,z(3)) +
(1 .-ycorr)*T(x(3),(YTC+1),z(3))

Tout(COUNT,4) = (xcorr*ycorr*T(XTC,YTC,z(4))+(1.-xcorr)*
ycorr*T((XTC+1 ),YTC,z(4))+xcorr*( 1 .-ycorr)*

143

00000

+ T(XTC,(YTC+1),z(4))+(1.-xcorr)*(1.-ycorr)*
+ T((XTC+1),(YTC+1),z(4)))

Tout(COUNT,5) = T(x(5), y(5), 2(5))

Tout(COUNT,6) = xcorr*T(XTC,y(6),z(6)) +
+ (1.-xcorr)*T((XTC+l),y(6),z(6))

Tout(COUNT,7) = ycorr*T(x(7),YTC,z(7)) +
+ (1 .-ycorr)*T(x(7),(YTC+1),z(7))

COUNT = COUNT+1
endif

close time loop

enddo

change parameters to account for larger timestep
fxx = axx*dt/(dx*dx)

fyy = ayy*dt/(dy*dy)

fzz = azz*dt/(dz*dz)

fxy = axy*dt/(2.*dx*dy)

hnode = (2.*hout*dt)/dz
backhnode = (2.*bkhout*dt)/dz
loop over rest of time

do 1 = 1, FINAL

here we compute the temperature field in the body

four top comers

Tnew(l,1,l) = fxx*(2.*T(2,1,1)-2.*T(1,1,l))+fyy*(2.*T(1,2,l)-
+ 2.*T(1 ,1,1))+fzz*(2.*T(1,1 ,2)-2.*T(1,1,1))+T(1,1,1)-hnode*T(1,1 ,1 )

Tnew(NX,1,1) = fxx*(2.*T(NX-1,1,1)-2.*T(NX,l,l))+fyy*(2.*
+ T(NX,2,])-2.*T(NX,1,1))+fzz*(2.*T(NX,l,2)-2.*T(NX,1,1))+T(NX,1,l )
+ -hnode*T(NX,1,1)

Tnew(l ,NY,1) = fxx*(2.*T(2,NY,l)-2.*T(1,NY,1))+fyy*(2.*
+ T(1,NY-l,1)-2.*T(1,NY,1))+fzz*(2.*T(1,NY,2)-2.*T(l,NY,1))+

144

 

 

 

+ T( 1 ,N Y, 1 )-hnode*T( 1 ,N Y, 1 )

Tnew(NX,NY,1) = fxx*(2.*T(NX-1 ,NY, l)-2.*T(NX,NY,1 ))+fyy*(2.*

T(NX,NY-1,1)-2.*T(NX,NY, l))+fzz*(2.*T(NX,NY,2)-2.*T(NX,NY,1))+

+
+ T(NX,NY, 1 )-hnode*T(NX,NY, 1)

two top edges

do i = 2, NX-l
Tnew(i,1,1)= fxx*(T(i+1,1,1)-2.*T(i,l ,1)+T(i-1,1,1))+
+ fyy*(2.*T(i,2,])-2.*T(i,1,1))+fzz*(2.*T(i,1,2)-2.*T(i,l,l))+

+ T(i, 1 ,1)-hnode*T(i,1,1)

Tnew(i,NY, l) = fxx*(T(i+1,NY,l)-2.*T(i,NY,l)+T(i-1,NY,1))+
+ fyy*(2.*T(i,NY-1,1)-2.*T(i,NY,1))+fzz*(2.*T(i,NY,2)-
+ 2.*T(i,N Y, 1 ))+T(i,N Y, 1 )-hnode*T(i,N Y, 1 )
enddo

two top edges

do j = 2, NY-l
Tnew(l,j,1) = fxx*(2.*T(2,j,1)-2.*T(1,j,l))+
+ fyy*(T(l,j+1,l)-2.*T(1,j,l)+T(1,j-l,1))+fzz*
+ (2.*T(l,j,2)—2.*T(l,j,1))+T(1,j,l)-hnode*T(1,j,1)

Tnew(NX,j, 1) = fxx*(2.*T(NX- l ,j,l )-2.*T(NX,j, 1))+
+ fyy*(T(NX,j+l ,1)-2.*T(NX,j,1)+T(NX,j-l,1))+fzz*

+ (2.*T(NX,j,2)-2.*T(NX,j,1))+T(NX,j,1)-hnode*T(NX,j,1)
enddo

top face
do i = 2, NX-l

do j = 2, NY-l

Tnew(i,j,l) = fxx*(T(i+1,j,1)-2.*T(i,j,1)+T(i-l,j,1))+

+ fyy*(T(i,j+1,1)-2.*T(i,j,l)+T(i,j-l,1))+fzz*(2.*T(i,j,2)-
+ 2.*T(i,j,1))+fxy*(T(i+1,j+l,l)-T(i+1,j-1,1)-T(i-1,j+l,l)+
+ T(i-l,j-l,1))+T(i,j,l)-hnode*T(i,j,l)

enddo
enddo

two side faces

doi=2,NX-1
do k = 2, NZ-l

145

Tnew(i, 1 ,k) = fxx*(T(i+1,1,k)-2.*T(i,1,k)+T(i-1,l,k))+

+ fyy*(2.*T(i,2,k)-2.*T(i,l,k))+fzz*(T(i,l,k+1)-2.*T(i,1,k)+
+ T(i,],k-1))+T(i,1,k)
Tnew(i,NY,k) = fxx*(T(i+1,NY,k)-2.*T(i,NY,k)+T(i-l,NY,k))+
+ fyy*(2.*T(i,NY-1 ,k)-2.*T(i,NY,k))+fzz*(T(i,NY,k+1)-2.*
+ T(i,NY,k)+T(i,NY,k-1))+T(i,NY,k)
enddo
enddo

two side faces

 

do j = 2, NY-l
do k = 2, NZ-l
Tnew( 1 ,j,k) = fxx*(2.*T(2,j,k)-2.*T(1 ,j,k))+fyy*
+ (T(l,j+l,k)-2.*T(l,j,k)+T(],j-l,k))+fzz*('1‘(l,j,k+1 )-2.*
+ T(1,j,k)+T(l,j,k-1))+T(1,j,k)
Tnew(NX,j,k) = fxx*(2.*T(NX- l ,j,k)-2.*T(NX,j,k))+fyy*
+ (T(NX,j+l,k)-2.*T(NX,j,k)+T(NX,j-l,k))+fzz*(T(NX,j,k+1)-
+ 2.*T(NX,j,k)+T(NX,j,k-l))+T(NX,j,k)
enddo
enddo

four side edges

do k = 2, NZ]

Tnew(1,1,k) = fxx*(2.*T(2,1,k)-2.*T(1,l,k))+fyy*
+ (2.*T(1,2,k)-2.*T(l,1,k))+fzz*(T(1,1,k+l)-2.*
+ T(l,1,k)+T(1,1,k-l))+T(l,1,k)

Tnew(N X, 1 ,k) = fxx*(2.*T(NX-1,1,k)-2.*T(NX,1 ,k))+fyy*
+ (2.*T(NX,2,k)-2.*T(NX,1,k))+fzz*(T(NX,1,k+1)-
+ 2.*T(NX,1,k)+T(NX,1,k-l))+T(NX,1,k)

Tnew(l,NY,k) = fxx*(2.*T(2,NY,k)-2.*T(l ,NY,k))+
+ fyy*(2.*T(1,NY-1,k)-2.*T(1,NY,k))+fzz*(T(1,NY,k+1)-2.*
+ T(l,NY,k)+T(l,NY,k-1))+T(1,NY,k)

Tnew(NX,NY,k) = fxx*(2.*T(NX-1,NY,k)-2.*T(NX,NY,k))+
+ fyy*(2.*T(NX,NY-l,k)-2.*T(NX,NY,k))+fzz*(T(NX,NY,k+1)-2.*
+ T(NX,NY,k)+T(NX,NY,k-l))+T(NX,NY,k)
enddo

model interior

146

 

do i = 2, NX-l

do j = 2, NY-l
do k = 2, NZ-l
Tnew(i,j,k) = fxx*(T(i+1,j,k)-2.*T(i,j,k)+T(i-l,j,k))+
+ fyy*(T(i,j+ 1 ,k)-2.*T(i,j,k)+T(i,j- l ,k))+fzz*
+ (T(i,j,k+1)-2.*T(i,j,k)+T(i,j,k-1))+fxy*(T(i+1,j+l ,k)-
+ T(i+1,j-1,k)-T(i—1,j+l,k)+T(i-1,j-1,k))+T(i,j,k)
enddo
enddo
enddo

bottom face

do i = 2, NX-l
do j = 2, NY-l
Tnew(i,j,NZ) = fxx*(T(i+1,j,NZ)-2.*T(i,j,NZ)+T(i-l,j,NZ))+

+ fyy*(T(i,j+],NZ)-2.*T(i,j,NZ)+T(i,j-1,NZ))+fzz*(2.*
+ T(i,j,NZ-1 )-2.*T(i,j,NZ))+fxy*(T(i+l ,j+1 ,NZ)-
+ T(i+1,j-l ,NZ)-T(i-1,j+l,NZ)+T(i-l,j-1,NZ))+T(i,j,NZ)-
+ backhnode*T(i,j,NZ)
enddo
enddo

two bottom edges

do i = 2, NX-l

Tnew(i,1,NZ) = fxx*(T(i+1,1,NZ)-2.*T(i,],NZ)+T(i-1,1,NZ))+
+ fyy*(2.*T(i,2,NZ)-2.*T(i,l ,NZ))+fzz*(2.*T(i,1,NZ-1)-
+ 2.*T(i,],NZ))+T(i,1,NZ)-backhnode*T(i,l,NZ)

Tnew(i,NY,N Z) = fxx*(T(i+1,NY,NZ)-2.*T(i,NY,NZ)+T(i-l,NY,NZ))+
+ fyy*(2.*T(i,NY-1,NZ)-2.*T(i,NY,NZ))+fzz*(2.*T(i,NY,NZ-1)-
+ 2.*T(i,NY,NZ))+T(i,NY,NZ)-backhnode*T(i,NY,NZ)
enddo

two bottom edges

do j=2, NY-l
Tnew(l,j,NZ) = fxx*(2.*T(2,j,NZ)-2.*T(l,j,NZ))+
+ fyy*(T(l,j+1,NZ)-2.*T(1,j,NZ)+T(],j-1,NZ))+fzz*
+ (2.*T(],j,NZ-l)-2.*T(1,j,NZ))+T(l,j,NZ)-backhnode*T(l,j,NZ)

Tnew(NX,j,NZ) = fxx*(2.*T(NX-1,j,NZ)-2.*T(NX,j,NZ))+

+ fyy*(T(NX,j-I-1,NZ)-2.*T(NX,j,NZ)+T(NX,j-l ,NZ))+fzz*(2.*
+ T(NX,j,N Z- 1 )-2.*T(NX,j,NZ))+T(NX,j,NZ)—backhnode*T(NX,j ,N Z)

147

enddo

four bottom comers

+
+
+

Tnew(1,l ,NZ) = fxx*(2.*T(2,l,NZ)-2.*T(l,1,NZ))+fyy*(2.*T(],2,NZ)—
2.*T(l,l,NZ))+fzz*(2.*T(1,1,NZ—1)-2.*T(1,1,NZ))+T(1,1,NZ)—
backhnode*T( l ,1 ,NZ)

Tnew(NX, 1 ,NZ) = fxx*(2.*T(NX-1,1,NZ)-2.*T(NX,1,NZ))+fyy*(2.*
T(NX,2,NZ)—2.*T(NX,1,NZ))+fzz*(2.*T(NX,1,NZ-l)-2.*T(NX,1,NZ))+
T(NX,1,NZ)-backhnode*T(NX,1 ,NZ)

Tnew( 1 ,NY,NZ) = fxx*(2.*T(2,NY,NZ)-2.*T(1,NY,NZ))+fyy*(2.*
T(l,NY-1,NZ)-2.*T(l,NY,NZ))+fzz*(2.*T(1,NY,NZ-1)-
2.*T(1,NY,NZ))+T(1,NY,NZ)-backhnode*T(1,NY,NZ)

Tnew(NX,NY,NZ) = fxx*(2.*T(NX-1,NY,NZ)-
2.*T(NX,NY,NZ))+fyy*(2.*T(NX,NY-1,NZ)-
2.*T(NX,NY,NZ))+fzz*(2.*T(NX,NY,NZ-1)-
2.*T(NX,NY,NZ))+T(NX,NY,NZ)-backhnode*T(NX,NY,NZ)

send Tnew to T once calculations are complete

T = Tnew

if the calculated time equals the desired output time, write the
temperature at the points of interest to a separate array

time = dt*l + 2.*flashtime
outtime = ( l .lfliz)*COUNT

if((time+del).ge.outtime) then

+

timeout(COUNT) = time
Tout(COUNT,1) = T(x(1), y(l), 2(1))

Tout(COUNT,2) = xcorr*T(XTC,y(2),z(2)) +
(1 .-xcorr)*T((XTC+l),y(2),z(2))

Tout(COUNT,3) = ycorr*T(x(3),YTC,z(3)) +
(1 .-ycorr)*T(x(3),(YTC+l ),z(3))

Tout(COUNT,4) = (xcorr*ycorr*T(XTC,YTC,z(4))+(1 .-xcorr)*
ycorr*T((XTC+1),YTC,z(4))+xcorr*(l .-ycorr)*
T(XTC,(YTC+1),z(4))+(1.-xcorr)*(l.-ycorr)*
T((XTC+1),(YTC+1),z(4)))

148

00000

0

Tout(COUNT,5) = T(x(5), y(5), 2(5))

Tout(COUNT,6) = xcorr*T(XTC,y(6),z(6)) +
+ (1 .-xcorr)*T((XTC+l),y(6),z(6))

Tout(COUNT,7) = ycorr*T(x(7),YTC,z(7)) +
+ (1 .-ycorr)*T(x(7),(YTC+1),z(7))

outtime = outtime + (1 ./f1‘lZ)
COUNT = COUNT+1
endif

close time loop

enddo

nreadm = COUNT-l

handle program output

send output temperature history to the MO program

do i = 1, nreadm

j = i+nreadm

k = i+2*nreadm

l = i+3*nreadm
set to center TC

tmod(i) = Tout(i,1)
set to old X TC

tmod(i) = Tout(i,2)
set to old Y TC

tmod(k) = Tout(i,3)
set to off axis sensor

tmod(l) = Tout(i,4)
enddo

if we are running the forward problem, add noise to simulate
experimental data.

if (noiseswitch.eq. 1) then
call SYSTEM_CLOCK(idum)
if (idum.eq.0) pause 'idum = 0; restart'
write(*,*) 'idum = ',idum
do i = 1,nreadm
do j = 1,5

149

0 ()(5 0

95

Tout(i,j) = Tout(i,j) + noise*gasdev(idum)
enddo

enddo
elseif(noiseswitch.eq.0) then

confinue
else

pause 'error with noiseswitch'
endif

open output file and write out temperature history at the points
of interest

open (unit=92, file="modelout.txt", status = "unknown")

do 1 = 1, nreadm
write(92,95) timeout(l),Tout(1, 1 ),Tout(l,2),Tout(l,3),

+ Tout(l,4),Tout(l,5)
enddo
format (1X,F7.3,2X,5F9.3)

close files, return to main program
close(92)

return
end

**************************

FUNCTION gasdev(idum)
integer idum
*uses ran2*
returns a normally distributed deviate with zero mean and unit
variance, using ran2(idum) as the source of uniform deviates
culled from Numerical Recipes in FORTRAN
integer iset
real(8) fac,gset,rsq,v 1 ,v2
save iset,gset
data iset/OI

if(iset.eq.0) then
v 1 =2.*ran2(idum)-1.
v2=2.*ran2(idum)-1.
rsq=v1**2+v2**2
if(rsq.ge. l ..or.rsq.eq.0.) goto 1

150

 

11

0

fac=sqrt(-2.*log(rsq)/rsq)
gset=vl *fac
gasdev=v2*fac
iset: 1

else
gasdev=gset
iset=0

endif

return
end

**************************

FUNCTION RAN2(IDUM)

from Numerical Recipes in FORTRAN

returns a uniform random deviate between 0. and l.
PARAMETER (M=7 14025,IA=1 366,IC=150889,RM=1 .4005 1 12E-6)
DIMENSION IR(97)
DATA IFF /0/
IF(IDUM.LT.0.0R.IFF.EQ.0)THEN

IFF=1

IDUM=MOD(IC-IDUM,M)

DO 11 J=l,97
IDUM=MOD(IA*IDUM+IC,M)
IR(J)=IDUM
CONTINUE

IDUM=MOD(IA*IDUM+IC,M)

IY=IDUM

ENDIF
J=1+(97*IY)/M
IF(J.GT.97.0R.J.LT.1)PAUSE
IY=IR(J)
RAN2=IY*RM
IDUM=MOD(IA*IDUM+IC,M)
IR(J)=IDUM
RETURN
END

****************************

end of subprograms, file

151

REFERENCES

152

REFERENCES

Aviles-Ramos, C., 2002, “Exact Solution of Heat Conduction in a Two-Domain
Parallelepiped with an Orthotropic Layer and its Application to the Estimation of
the Three-Dimensional Thermal Conductivity Tensor and Volumetric Heat
Capacity of the Orthotropic Layer,” Proceedings of the ASME IMECE 2002, New
Orleans, Louisiana, November 17-22 2002, IMECE2002-32424.

Beck J .V., and Arnold, K.J., 1977, Parameter Estimation in Engineering and Science,
John Wiley and Sons, New York.

Bennett, T.D., 2004, “Determining Anisotropic Film Thermal Properties Through
Harmonic Surface Heating With a Gaussian Laser Beam: A Theoretical
Consideration,” Journal of Heat Transfer 126, pp. 305-31 1.

Bradley, CB, 1943, “Measurement of thermal conductivity of thermal insulating
materials by means of guarded hot plate,” Refrigerating Engineering 45(5), pp.
337-341.

Broerman, A.W., Venerus, DC, and Schieber, J.D., 1999, “Evidence for the Stress-
Therrnal Rule in an Elastomer Subjected to Simple Elongation,” Journal of
Chemical Physics, 111(15), pp. 6965-6969.

Carslaw, HS, and J aeger, J.C., 1959, Conduction of Heat in Solids, 2"d ed., Oxford
University Press, Oxford, UK.

Chang, Y.P., and Tsou, R.C.H., 1977, “Heat Conduction in an Anisotropic Medium
Homogeneous in Cylindrical Regions - Unsteady State,” Journal of Heat
Transfer, 99, pp. 41-46.

Choy, C.L., Luk, W.H., and Chen, EC, 1978, “Thermal Conductivity of Highly Oriented
Polyethylene,” Polymer, 19(2), pp. 155-162.

Chu, F.I., Taylor, RE, and Donaldson, A.B., 1980, “Thermal Diffusivity Measurements
at High Temperatures by the Radial Flash Method,” Journal of Applied Physics,
51(1) pp. 336-341.

CRC Handbook of Chemistry and Physics, 71St ed., D. Lide, editor, CRC Press, Boston.

Davis, SE, and Wright, N.T., 2005, “An Optimal Experimental Design Approach for

Sensor Positioning in an Extended Flash Thermal Diffusivity Experiment”, in
preparation.

153

Degiovanni, A., 1977, “Diffusivite et methode flash (Thermal Diffusivity and the Flash
Method),” Revue Generale de Thermique, 16(185), pp. 420—441.

Donaldson, A.B., and Taylor, R.B., 1975, “Thermal Diffusivity Measurement by a Radial
Heat Flow Method,” Journal of Applied Physics, 46 (10), pp. 4584-4589.

Doss, D.J., 1998, “Thermal Diffusivity of Vascular Tissue and Polymers: Effects of
Finite Deformation and Direction,” Ph.D. dissertation, University of Maryland
Baltimore County, Baltimore, MD.

Doss, D.J., and Wright, N.T., 2000, “Simultaneous Measurement of the Orthogonal
Components of Thermal Diffusivity in PVC Sheet,” Journal of Heat Transfer,
122(1), pp. 27-32.

Douglas, J ., 1962, “Alternating Direction Methods for Three Space Variables,”
Numerische Mathematik, 41, pp. 41-63.

Dowding, K.J., and Blackwell, BE, 1999, “Design of Experiments to Estimate
Temperature Dependent Thermal Properties,” Proceedings of ASME 3”
Conference on Inverse Problems in Engineering: Theory and Practice, pp. 509-
5 18.

Emery, A.F., 2001, “Using the Concept of Information to Optimally Design Experiments
with Uncertain Parameters,” Journal of Heat Transfer, 123, pp. 593-600.

Emery, A.F., and Fadale, T.D., 1997, “The Effect of Imprecisions in Thermal Sensor
Location and Boundary Conditions on Optimal Sensor Location and Experimental
Accuracy,” Journal of Heat Transfer, 119, pp. 661-665.

Emery, A.F., Nenarokomov, A.V., and Fadale, T.D., 2000, “Uncertainties in Parameter
Estimation: The Optimal Experiment Design,” International Journal of Heat and
Mass Transfer, 43, pp. 3331-3339.

Fadale, T.D., Nenarokomov, A.V., and Emery, A.F., 1995, “Two Approaches to Optimal
Sensor Locations,” Journal of Heat Transfer, 117, pp. 373-379.

Graham, 8., McDowell, D.L., and Dinwiddie, R.B., 1999, “In-Plane Thermal Diffusivity
Measurements of Orthotropic Materials,” Thermal Conductivity 24, RS. Gaal, ed.,
Technomic, Lancaster, PA, pp. 241-252.

Graham, 8., McDowell, D.L., and Dinwiddie, R.B., 1999, “Multidimensional Flash

Diffusivity Measurements of Orthotropic Materials,” International Journal of
Thermophysics, 20(2), pp. 691-707.

154

Haji-Sheikh, A., and Beck, J.V., 2002, “Temperature Solution in Multi-Dimensional
Multi-Layer Bodies,” International Journal of Heat and Mass Transfer, 45, pp.
1865- 1877.

Hsieh, M., and Ma, Q, 2002, “Analytical Investigations for Heat Conduction Problems in
Anisotropic Thin-Layer Media With Embedded Heat Sources,” International
Journal of Heat and Mass Transfer, 45, pp. 4117-4132.

Huang, SC, and Chang, Y.P., 1984, “Anisotropic Heat Conduction With Mixed
Boundary Conditions,” Journal of Heat Transfer, 106, pp. 646-648.

Huang, X., 2004, “Comparison of Thermal Diffusivity Measurement Techniques,” M.S.
thesis, Michigan State University, East Lansing, MI.

Incropera, F., and DeWitt, D., 1990, Fundamentals of Heat Transfer, 3'“ edition, John
Wiley and Sons, New York.

Kato, H, Baba, T, and Okaji, M, 2001, “Anisotropic Thennal-Diffusivity Measurements
by a New Laser-Spot-Heating Technique,” Measurement Science and
Technology, 12, pp. 2074-2080.

LeGall, E.L., 2001, “Thermal Diffusivity of Elastomers Subject to Finite Biaxial
Deformation,” Ph.D. dissertation, University of Maryland Baltimore County,
Baltimore, MD.

LeGall, E.L., and Wright, N.T., 2000, “Changes in Thermal Diffusivity of Rubber Due to
Mechanical Preconditioning”, paper presented at the Fourteenth Symposium on
Thermophysical Properties, June 2000, Boulder Colorado, pp. 1-22.

Le Niliot, C., Rigollet, F., and Petit, D., 2000, “An Experimental Identification of Line
Heat Sources in a Diffusive System Using the Boundary Element Method,”
International Journal of Heat and Mass Transfer, 43, pp. 2205-2220.

Ma, C. and Chang, S., 2004, “Analytical Exact Solutions of Heat Conduction Problems
for Anisotropic Multi-Layered Media,” International Journal of Heat and Mass
Transfer, 47, pp. 1643-1655.

Maillet, D., Lachi, M. and Degiovanni, A., 1990, “Simultaneous Measurements of Axial
and Radial Thermal Diffusivities of an Anisotropic Solid in Thin Plate:
Application to Multi-Layered Materials,” Thermal Conductivity 21, Cremers, OJ.
and Fine, H.A., eds., Plenum Press, New York, pp. 91-107.

Mandelis, A., 1985, “Hamilton-Jacobi Formulation and Quantum Theory of Thermal

Wave Propagation in the Solid State,” Journal of Mathematical Physics, 26(10),
pp. 2676-2683.

155

Mandelis, A., Funak, F., and Munidasa, M., 1996, “Generalized Methodology for
Thermal Diffusivity Depth Profile Reconstruction in Semi-infinite and Finitely
Thick Inhomogeneous Solids,” Journal of Applied Physics, 80(10), pp. 5570-
5578.

Matthews, G., 1996, PVC: Production, Properties and Uses, The University Press,
Cambridge.

McMasters, R.L., 1997, “Modeling Heat Transfer for Parameter Estimation in Flash
Diffusivity Experiments,” Ph.D. dissertation, Michigan State University, East
Lansing, MI.

McMasters, R.L., Beck, J.V., Dinwiddie, R.B., and Wang, H., 1999, “Accounting for
Penetration of Laser Heating in Flash Thermal Diffusivity Experiments,” Journal
of Heat Transfer, 121, pp. 15-21.

Mejias, M.M., Orlande, H.R.B., and Ozisik, M.N., 1999, “Design of Optimum
Experiments for the Estimation of the Thermal Conductivity Components of
Orthotropic Solids,” Hybrid Methods in Engineering, 1, pp. 37-53.

Mejias, M.M., Orlande, H.R.B., and Ozisik, M.N., 2000, “On the Choice of Boundary
Conditions for the Estimation of the Thermal Conductivity Components of
Orthotropic Solids,” Proceedings of NH TC ’00, 34m National Heat Transfer
Conference, Pittsburgh, Pennsylvania, August 20-22, 2000, NHTC2000-12062.

Mejias, M.M., Orlande, H.R.B., and Ozisik, M.N., 2003, “Effects of the Heating Process
and Body Dimensions on the Estimation of the Thermal Conductivity

Components of Orthotropic Solids,” Inverse Problems in Engineering, 11(1), pp.
75-89.

Milosevic, ND, and Raynaud, M., 2004, “Analytical Solution of Transient Heat
Conduction in a Two-Layer Anisotropic Cylindrical Slab Excited Superficially by
a Short Laser Pulse,” International Journal of Heat and Mass Transfer, 47 , pp.
1627-1641.

Mulholland, GP, and Gupta, BR, 1977, “Heat Transfer in a Three-Dimensional
Anisotropic Solid of Arbitrary Shape,” Journal of Heat Transfer, 99, pp. 135-137.

Mzali, F., Sassi, L., Jemni, A., Ben Nasrallah, S., and Petit, D., 2004, “Optimal
Experiment Design for the Identification of Thermo-Physical Properties of

Orthotropic Solids,” Inverse Problems in Science and Engineering, 12(2), pp.
193-209.

Onsager, L., 1931, “Reciprocal Relations in Irreversible Processes. 1,” Physical Review,
37(4) pp. 405-426.

156

Ozisik, M.N., 1993, Heat Conduction, 2'“! ed., John Wiley and Sons, Inc., New York.

Parker, W.J., Jenkins, R.J., Butler, CR, and Abbott, G.L., 1961, “Flash Method of
Determining Thermal Diffusivity, Heat Capacity, and Thermal Conductivity,”
Journal of Applied Physics, 32(9), pp. 1679-1684.

Petit, D., and Hachette, R., 1998, “Model Reduction in Linear Heat Conduction: Use of

Interface Fluxes for the Numerical Coupling,” International Journal of Heat and
Mass Transfer, 41, pp. 3177-3189.

Philippi, 1., Batsale J .C., Maillet, D., and Degiovanni, A., 1995, “Measurement of
Thermal Diffusivities Through Processing of Infrared Images,” Review of
Scientific Instruments, 66(1), pp. 182-192.

Poon, K.C., Tsou, R.C.H., and Chang, Y.P., 1979, “Solution of Anisotropic Problems of
First-Class by Coordinate Transformation,” Journal of Heat Transfer, 101, pp.
340-345.

Powers, J.M., 2004, “On the Necessity of Positive Semi-Definite Conductivity and
Onsager Reciprocity in Modeling Heat Conduction in Anisotropic Media,”
Journal of Heat Transfer, 126, pp. 670-675.

Pozrikidis, C., 1998, Numerical Computation in Science and Engineering, New York:
Oxford University Press, pp. 526-550.

Press, W.H., Teukolsky, S.A., Vetterling, WT, and Flannery, BR, 1992, Numerical
Recipes in FORTRAN: The Art of Scientific Computing, 2"d ed., Cambridge
University Press, New York.

Pron, H., and Bissieux, C., 2004, “3-D Thermal Modelling Applied to Stress-Induced
Anisotropy of Thermal Conductivity,” International Journal of Thermal Sciences,
43, pp. 1161-1169.

Sawaf, B., and Ozisik, M.N., 1995, “Determining the Constant Thermal Conductivities of
Orthotropic Materials by Inverse Analysis,” International Communications on
Heat and Mass Transfer, 22(2), pp. 201-211.

Taktak, R., 1992, “Design and Validation of Optimal Experiments for Estimating
Thermal Properties of Composite Materials,” Ph.D. dissertation, Michigan State
University, East Lansing, MI.

Taktak, R., Beck, J.V., and Scott, ER, 1993, “Optimal Experimental Design for

Estimating Thermal Properties of Composite Materials,” International Journal of
Heat and Mass Transfer, 36(12), pp. 2977-2986.

157

Torres, J.H., Springer, T.A., Welch, A.J., and Pearce, J .A., 1990, “Limitations of a
Thermal Camera in Measuring Surface Temperature of Laser Irradiated Tissues,”
Lasers in Surgery and Medicine, 10, pp. 510-523.

Venerus, D.C., Schieber, J.D., Iddir, H., Guzman, J .D., and Broerman, A.W., 1999,
“Relaxation of Anisotropic Thermal Diffusivity in a Polymer Melt Following Step
Shear Strain,” Physical Review Letters, 82(2), pp. 366-369.

Venerus, D.C., Schieber, J .D., Iddir, H., Guzman, J.D., and Broerman, A.W., 1999,
“Measurement of Thermal Diffusivity in Polymer Melts Using Forced Rayleigh

Light Scattering,” Journal of Polymer Science Part B: Polymer Physics, 37, pp.
1069- 1078.

Venerus, D.C., Schieber, J.D., Iddir, H., Guzman, J ., and Broerman, A., 2001,
“Anisotropic Thermal Diffusivity Measurements in Defonning Polymers and the

Stress-Thermal Rule,” International Journal of Thermophysics, 22(4) pp. 1215-
1225.

Videcoq, E., and Petit, D., 2001, “Model Reduction for the Resolution of

Multidimensional Inverse Heat Conduction Problems,” International Journal of
Heat and Mass Transfer, 44, pp. 1899-1911.

158