LIBRARY

Michigan State
Unlvarllty

 

 

 

PLACE ll RETURN BOX to romovo this checkout from your mood.
TO AVOID FINES Mum on or baton date duo.

DA'I; DUE DATE DUE DATE DUE ll

 

 

 

NEW-r69 ]
: l

 

 

 

 

 

 

 

 

 

 

———l
A;
:1

 

 

1
l

MSU I. An Affirmative Action/Equal Opponunuy Institution
Wanna-9.1

 

MULTI-DIMENSIONAL ESTIMATION OF THERMAL PROPERTIES
AND SURFACE HEAT FLUX USING EXPERIMENTAL DATA AND A
SEQUENTIAL GRADIENT METHOD
By

Kevin J Dowding

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

1997

ABSTRACT

MULTI-DIMENSIONAL ESTIMATION OF THERMAL PROPERTIES
AND SURFACE HEAT FLUX USING EXPERIMENTAL DATA AND A
SEQUENTIAL GRADIENT METHOD

By

Kevin J Dowding

Inverse theory is an emerging field of study with application to a diverse range of
problems. Inverse thermal problems are the focus of this dissertation; specifically, the
parameter estimation problem and inverse heat conduction problem (H-ICP) are investi-
gated. Although one-dimensional inverse thermal problems have been widely investigated,
multi-dimensional problems are beginning to receive an increasing amount of attention.
One- and two-dimensional cases are addressed for both the noted inverse problems,
including an experimental application.

Parameter estimation techniques are applied to estimate the thermal properties of a
carbon-carbon composite from transient experiments. Properties are determined as a func-
tion of temperature and direction relative to the fiber orientation. The thermal conductivity
is assumed to be orthotropic, varying in the direction normal and parallel to the fibers; the
volumetric heat capacity is assumed isotropic. Thermal properties from room temperature
up to 500C are obtained. The thermal conductivity normal to the fiber is found to be less
than one-tenth of the thermal conductivity parallel to the fiber. Agreement within 7% is

demonstrate between independent one- and two-dimensional results.

A sequential-in-time implementation is proposed for a conjugate gradient method,
utilizing an adjoint equation approach, to solve the IHCP. Because the IHCP is generally
ill-posed, Tikhonov regularization is included to stabilize the solution. The proposed
sequential method benefits from the efficiency and on-line capabilities of a sequential
implementation, without requiring a priori information about the (unknown) surface heat
flux. Aspects of the sequential gradient method are discussed and examined. Several
promising features of the sequential gradient method are noted. Simulated one- and two-
dimensional test cases are presented to study the sequential implementation. Numerical
solutions are obtained using a finite difference procedure. Results indicate the sequential
implementation has accuracy comparable to a standard whole domain solution, but in cer-
tain cases requires significantly more computational time. Methods to improve the compu-

tational requirements, which make the method competitive, are presented.

to my wife

iv

ACKNOWLEDGMENTS

The input and help of the my Ph.D. guidance committee, James Beck, Patti Lamm,
Alex Diaz, and Craig Somerton, is acknowledged. Special thanks are due Patti Lamm for
her patience working with an engineer, help with understanding the gradient methods, and
availability to discuss the research. To my advisor James Beck, I thank for his encourage-
ment and introduction to the field of inverse problems. His enthusiasm and dedication
were inspiring, not only in helping realize this work, but as a teacher and researcher.

I have benefited from interaction with several fellow students and researchers. Input
and help from Arafa Osman, Bob McMasters, Heidi Relyea, and Matt White is appreci-
ated.

Completing this work would not have been possible if it were not for the endless sup-

port my parents. To my wife I’m indebted for her understanding, support, and patience.

TABLE OF CONTENTS

 

 

 

 

 

 

LIST OF TABLES ix
LIST OF FIGURES Y
NOMENCLATURE xiv
Greek ................................................................................................................................ xvi
Subscripts ........................................................................................................................ xvii
Superscripts .................................................................................................................... xviii
Chapter 1
INTRODUCTION 1
Chapter 2
CARBON -CARBON THERMAL PROPERTIES: ONE-DIMENSIONAL EXPERI-
MENTS 6
2-1.0 Introduction .............................................................................................................. 6
2-l.l Problem Description ............................................................................................................ 8
2-1.2 Literature review .................................................................................................................. 9
2-2.0 Experimental Aspects ............................................................................................ 11
2-2.1 Experimental Description .................................................................................................. 11
2-2.2 Experimental Design .......................................................................................................... 19
2-3.0 Analysis Procedures ............................................................................................... 21
2-4.0 Results and Discussion .......................................................................................... 23
2-4.1 Effective Properties for Mica Heater and Insulation .......................................................... 24
2-4.2 Effective Properties of Carbon-Carbon .............................................................................. 27
2-5.0 Experimental Uncertainty ...................................................................................... 36
Chapter 3
CARBON-CARBON THERMAL PROPERTIES: TWO-DIMENSIONAL EXPERI-
MENTS 38
3-1.0 Introduction ............................................................................................................ 38
3-1.1 Motivation .......................................................................................................................... 38
3-2.0 Experimental Aspects ............................................................................................ 40

vi

3-3.0 Analysis Procedures ............................................................................................... 43

 

 

3-4.0 Results and Discussion .......................................................................................... 44
3-4.1 Experimental Data ............................................................................................................. 45
3-4.2 Estimated Thermal Properties ............................................................................................ 47
3-4.3 Comparison of One- and Two-Dimensional Results ......................................................... 60
3-5.0 Experimental Uncertainty ...................................................................................... 63
Chapter 4
SOLUTION OF IHCP USING A GRADIENT METHOD WITH ADJOINT EQUA-
TION APPROACH 65
4-1.0 Introduction ............................................................................................................ 65
4-2.0 Problem Statement of the Multi-Dimensional IHCP ............................................. 70
4-3.0 Gradient Calculation .............................................................................................. 73
4-4.0 Sensitivity Equations ............................................................................................. 75
4-5.0 Adjoint Equations .................................................................................................. 77
4-6.0 Optimization Method ............................................................................................. 84
4-7.0 Implementation of the Gradient Method ................................................................ 87
4-7.1 Whole Domain ................................................................................................................... 87
4-7.2 Sequential Implementation ................................................................................................ 91
4-8.0 Parameterizing the Solution (Finite-Dimensional Problem) .................................. 96
4-9.0 Summary ................................................................................................................ 99
Chapter 5
APPLICATION OF THE SEQUENTIAL GRADIENT METHOD: ONE-DIMEN-
SIONAL IHCP 101
5-1.0 Introduction .......................................................................................................... 101
5-2.0 Advantages of a Sequential Method and Numerical Solution of the IHCP ......... 103
52.] Numerical Solution .......................................................................................................... 103
5-2.2 Sequential Solution of IHCP ............................................................................................ 104
5-3.0 Simulated Test Cases ........................................................................................... 107
5-4.0 1D Results - Simulated Measurements ................................................................ 115
5-4.1 Exact Data (no Measurement Errors) .............................................................................. l 15
5-4.2 Corrupted Data (Measurement Errors) ............................................................................ 126
5-4.3 Results - Experimental Measurements ............................................................................. 140
Chapter 6
APPLICATION OF THE SEQUENTIAL GRADIENT METHOD: TWO-DIMEN-
SIONAL IHCP 145

 

vii

 

 

 

 

6-1.0 Introduction .......................................................................................................... 145
6-2.0 Simulated Temperature Data ................................................................................ 146
6-3.0 2D Results - Simulated Measurements ................................................................ 149
6-3.1 Exact Data (No Measurement Errors) .............................................................................. 149
6-3.2 Corrupted Data (Measurement Errors) ............................................................................ 159
6-4.0 Results - Experimental Measurements ................................................................. 175
Chapter 7
SUMMARY AND CONCLUSIONS 185
Chapter 8
RECOMMENDED FUTURE WORK 187
APPENDIX A
FINITE CONTROL VOLUME METHOD 189
A-1.0 Introduction .......................................................................................................... 189
A-2.0 Alternating Direction Implicit Method (ADI) ..................................................... 195
A-2.l ADI Equations (n+l/2) .................................................................................................... 196
A-2.2 ADI Equations (n+1) ....................................................................................................... 198
A-2.3 Summary of ADI Equations ............................................................................................. 200
LIST OF REFERENCE 201

viii

Table 2-1

Table 2-2

Table 2-3

Table 2-4

Table 2-5

Table 3-1

Table 3—2

Table 5-1

Table 5-2

Table 5-3

Table 5-4

Table 6-1

Table 6-2

Table 6-3

Table 6-4

LIST OF TABLES

Sensor locations in experimental apparatus ............................................... 14
Thermal properties of the ceramic insulation ............................................ 25
Effective thermal properties for the mica heater assembly ........................ 25

Thermal properties estimated for the carbon-carbon composite from one-
dimensional experiments and analysis ....................................................... 27

Experimental uncertainty for the two-dimensional experiments ............... 37

Thermal properties estimated for the carbon-carbon composite from two-

dimensional experiments and analysis ....................................................... 48
Experimental uncertainty for the two-dimensional experiments ............... 63
Estimation results for exact simulated data ............................................. 117
Estimation results for simulated data corrupted with errors .................... 130

Estimation results for prescribing a constant heat flux over the sequential
interval ..................................................................................................... 138

Estimation results for analysis of experimental case ............................... 144

Estimation results for the two-dimensional IHCP with exact simulated
data ........................................................................................................... 155

Estimation results for two-dimensional IHCP with simulated data
corrupted with random errors ................................................................... 162

Estimation results for two-dimensional IHCP with simulated data
corrupted with random errors using prior information ............................ 174

Estimation results for two-dimensional IHCP with experimentally
measured heat flux ................................................................................... 183

ix

Figure 2-1

Figure 2-2

Figure 2-3

Figure 2-4
Figure 2-5
Figure 2-6
Figure 2-7
Figure 3-1

Figure 3-2

Figure 3-3

Figure 3-4

Figure 3-5a

Figure 3-5b

Figure 3-5c

Figure 3-6a

LIST OF FIGURES

Schematic of experimental set-up used for estimating thermal properties
of carbon-carbon composite. For one-dimensional experiments all heaters
are activated. .............................................................................................. 12

Heat transfer model for one-dimensional experiment ............................... 15

Temperature dependence of thermal properties estimated from one-
dimensional experiments ........................................................................... 29

Experimental data for one-dimensional case 1023#508.2 ......................... 30

Sequential parameter estimates for one-dimensional case 1023#508.2 ...32

Temperature residuals for one-dimensional case 1023#508.2 ................... 33
Sensitivity coefficients for one-dimensional case 1023#508.2 .................. 35
Heat transfer model for two-dimensional experiments .............................. 41

Experimental data for two-dimensional case 1022$297, see Table 2-1 for
sensor locations .......................................................................................... 45

Temperature dependence of estimated thermal properties from two-
dimensional experiments ........................................................................... 49

Sequential parameter estimates for two-dimensional case 1022$297 ....... 50

Temperature residuals for two-dimensional case 1022$297 on surface
below active heater ..................................................................................... 52

Temperature residuals for two-dimensional case 1022$297 on the heated
surface, but not on the are of the active heater ........................................... 53

Temperature residuals for two-dimensional case 1022$297 at the
specimen/insulation interface ..................................................................... 54

Normalized sensitivity coefficient for volumetric heat capacity for two-
dimensional case 1022$297 ....................................................................... 57

Figure 3-6b

Figure 3-6c

Figure 3-7

Figure 4-1

Figure 5-1

Figure 5-2

Figure 5-3

Figure 5-4

Figure 5-5a

Figure 5-5b

Figure 5-5c

Figure 5-6

Figure 5-7

Figure 5-8

Figure 5-9

Figure 5-10

Normalized sensitivity coefficient for thermal conductivity in y-direction
for two-dimensional case 102233297 .......................................................... 58

Normalized sensitivity coefficient for thermal conductivity in x-direction
for two-dimensional case 1022$297 .......................................................... 59

Comparison of thermal properties estimated from one— and two-
dimensional experiments a) thermal conductivity and b) volumetric
specific heat ................................................................................................ 62

Schematic of multi-dimensional general IHCP ......................................... 70

Simulated temperature and heat flux data for the one-dimensional
triangular heat flux ................................................................................... 112

Simulated temperature and heat flux data for the one-dimensional step
heat flux test case ..................................................................................... 1 13

Simulated temperature and heat flux data for the one-dimensional
sinusoidal heat flux test case .................................................................... 1 14

Estimated heat flux for one-dimensional triangular test case with exact
data ........................................................................................................... 116

Estimated one-dimensional triangular heat flux with no measurement
errors for sequential interval beginning at time 0.12 for r = 5 ................. 123

Estimated one-dimensional triangular heat flux with no measurement
errors for sequential interval beginning at time 0.12 for r = 8 ................. 124

Estimated one-dimensional triangular heat flux with no measurement
errors for sequential interval beginning at time 0.12 for r = 10 ............... 125

Estimated heat flux for one-dimensional triangular test case with
measurement errors .................................................................................. 127

Estimated heat flux for one-dimensional step test case with measurement
errors ........................................................................................................ 128

Estimated heat flux for one-dimensional sinusoidal test case with
measurement errors .................................................................................. 129

Procedure for specifying the initial condition using converged estimates
from the previous sequential interval ....................................................... 135

Estimated heat flux for one-dimensional triangular case with the heat flux
held constant over future sequential interval ........................................... 139

xi

Figure 5-11
Figure 6-1

Figure 6-2

Figure 6—3

Figure 6-4a

Figure 6-4b

Figure 6-5a

Figure 6-5b

Figure 6-6a

Figure 6—6b

Figure 6-7

Figure 6-8

Figure 6—9a

Figure 6-9b

Figure 6-9c

Figure 6-9d

Estimate heat flux for one-dimensional experimental case 1010#30.1 ..143
Two-dimensional geometry for simulated test cases ............................... 147
Simulated temperature and surface heat flux for a two-dimensional

triangular heat flux ................................................................................... 150

Simulated temperature and surface heat flux for a two-dimensional step
heat flux .................................................................................................... 151

Estimated surface two-dimensional heat flux for triangular test case using
data without measurement errors. Whole domain solution ..................... 153

Estimated surface two-dimensional heat flux for triangular test case using
data without measurement errors. Sequential solution with 1:6 ............. 154

Estimated surface heat flux for triangular test case using data corrupted
with measurement errors (0’ = 0.0018°C ). Whole domain solution. 160

Estimated surface heat flux for triangular test case using data corrupted
with measurement errors (6 = 0.0018°C ). Sequential solution r=6 161

Estimated surface heat flux for step test case using data corrupted with
measurement errors (0’ = 0.0025°C ). Whole domain solution ............. 164

Estimated surface heat flux for step test case using data corrupted with
measurement errors (0' = 0.0025°C ). Sequential solution r=15 ........... 165

Estimated surface heat flux for triangular test case using data corrupted
with measurement errors ( 0' = 0.0018°C ). Sequential solution r=15 168

Estimated surface heat flux for triangular test case using data corrupted

with measurement errors (6 = 0.0018°C ). Sequential solution with $6
using prior information. ........................................................................... 173

Measured surface heat for experimental case to estimate thermal
properties of carbon-carbon ..................................................................... 177

Estimated surface heat flux for experimental case. Whole domain ......... 180

Estimated surface heat flux for experimental case. Sequential solution
with r-—10. ................................................................................................. 181

Estimated surface heat flux for experimental case. Combined function
specification regularization method with r=6. ......................................... 182

xii

Figure A-l Computational nodes and control surfaces .............................................. 191
Figure A-2 Typical energy balance for an interior node ............................................. 191

Figure A-3 Finite control volumes along the boundary of the domain ...................... 194

xiii

”f."

NOMENCLATURE

dimension of two-dimensional simulated case [m]
area [m2]

sub-diagonal of tri-diagonal matrix

dimension of two-dimensional simulated case [m]
parameter vector

estimated parameter vector

sup-diagonal of tri-diagonal matrix

specific heat [J/kgC]

main diagonal of tri-diagonal matrix

location for temperature sensor j [m]

temperature residual sensor j [°C]

function space

nonhomogeneous term for boundary surface [1“,]
convection coefficient [W/(m2°C)]

boundary coefficient

number of components retained in sequential implementation
current [Amps]

objective function [°C2 sec]

sum-of-squares term in objective function [°Czsec]

regularization (Tikhonov) term if objective function [°Czsec]

xiv

*0 23: =>

”6
a

qprr'or

qest

thermal conductivity [W/m° C]
boundary coefficient
length [m]
function space of all “square integrable” functions
number of temporal components for estimated heat flux on [F4]
outward pointing unit normal vector
number of measurement times
number of estimated parameters
number of parameters
search direction [—o—:——:|
(W/m )m
number of spatial components for heat flux on [F4]
heat flux [W/mz]
estimated heat flux [W/mz]
estimated heat flux with measurement errors [W/mz]
prior information for heat flux, equation (4—2.4) [W/mz]
stored flux for sequential solution [W/mz]
two-dimensional coordinate vector [m]
number of future time steps
sum-of-squares function [°C2]
temperature error, equation (5-4.2) [°C]
basis for finite-dimensional heat flux, equation (4-7.6)
time [sec]
heating time [sec]
temperature [°C]

XV

?<€<
><

?<l
><l

7i

7i, k

At

M

“i

calculated temperature [° C]

prior information vector

weighting matrix for prior information
velocity vector [m/s]

voltage [V]

 

. . . 1
weighting matrix for sensors [ 2]
°C
sensitivity coefficient, matrix
normalized sensitivity coefficient, matrix [°C]

measured temperature, vector [°C]

thermal diffusivity [mzls]

. . . °C 2 1
Trkhonov regularization parameter I: — —]

(W/mz) m

conjugate gradient search direction parameter, equation (4-6.10)
coefficient in tri-diagonal solution, equation (5-2.3b)
coefficient in tri-diagonal solution, equation (5-2.3d)
weighting for finite-dimensional heat flux, equation (4-7.5)
boundary surface i for domain (2
optimality criteria
time step [sec]
update to heat flux [W/mz]
convergence tolerance
temperature error [°C]

Lagrange multiplier, equation (4-4.3)

xvi

9 sensitivity function [°C]

 

 

0 adapted sensitivity function [(“jfnzfa

7t Lagrange multiplier, equation (443)

A Lagrangian function, equation (4-4.3)

p density [kg/m3]

fin step size, equation (4-6.l3) [(vK/gszm]

0’ standard deviation of temperature noise [°C]

60 estimated error in heat flux, no measurement errors, equation (5-4.3) [°C]
65‘ estimated error in heat flux, random measurement errors, equation (5-4.7) [°C]
6,, estimated error in heat flux, experimental data [°C]

67:. estimated error in measured temperature [°C]

1: relaxation parameter for Residual Principle

(1)1, (1)2, (b3 one-dimensional solutions for basic heat flux functions [°C]

 

4),- basis parameter on time, equation (4-7.1)
(pk basis function on space, equation (4-7. 1)
. . . . . °C 1
ad ornt function La ran e multr lrer [ —]
\V J ( g g p ) (W/m2)m
Q two-dimensional domain
Subscripts
0 initial
cc carbon-carbon composite material
f final
ins insulating material

kapton kapton heater assembly

xvii

mica mica heat assembly

x direction parallel to the fiber direction
y direction normal to the fiber direction
Superscripts

e effective

n, k iteration index

+ dimensionless variable

xviii

Chapter 1

INTRODUCTION

Study of direct (well-posed) problems has progressed for nearly two centuries. In heat
conduction, a direct problem represents a situation when the coefficients of the describing
differential equation, such as thermal conductivity, density, and specific heat, in addition
to the magnitude and location of internal energy generation (if it exist), are specified. The
geometry of the body, boundary conditions, and initial conditions are also specified. Then
given this specified information, the temperature at any point in the body is calculated.
Direct problems have been examined for nonlinear and multi-dimensional cases as well as
irregular geometries. Study of direct problems is widely described in the literature, includ-
ing issues of uniqueness and stability. Although analytical solution methods, Beck, et al.
(1992) and Ozisik (1993), are restricted to mainly linear problems with a regular geome-
try, many numerical methods are available for solving more complex direct problems,
Minkowycz et al. (1988).

Direct problems are typically in a class called well-posed problems because their
solutions exist, are unique, and depend continuously on the data (coefficients of describing
equation, internal energy generation, and boundary and initial conditions), Lamm (1993).
We can think of the direct problem as a cause-effect relationship. Given the cause (coeffi-

cients of describing equation, energy generation, and boundary and initial conditions), the

2

effect (temperature) can be determined. In contrast to the direct problem, an inverse prob-
lem attempts to determine the cause for a given efiect; the inverse of the direct problem.
That is, given the effect (temperature), the cause (coefficients of describing equation,
energy generation, or boundary and initial conditions) is determined. Inverse problems are
typically ill-posed because they do not satisfy the criteria for a well-posed problem. Prob-
lems that do not satisfy the criteria for being well-posed were initially thought to be of lit-
tle practical and physical significance. As it turns out, inverse problems are not only
important, but are applicable to a wide range of practical problems.

Inverse thermal methods have numerous applications. Specific examples include: the
determination of the outer surface heat flux during re-entry of a space vehicle, estimating
surface conditions at the exhaust of a rocket or jet engine, determining conditions in the
cylinder of an internal combustion engine, analysis of quenching in material processing,
industrial process control, optimal design of experiments, estimating thermophysical
properties of materials, analysis of casting processes, etc. Two popular applications of
inverse methods are concerned with the estimation of surface conditions from temperature
measurements, typically called the inverse heat conduction problem (IHCP) and the esti-
mation of thermophysical properties or parameter estimation (PE). Both of these prob-
lems are addressed in this dissertation.

The study of inverse problems is relatively recent, particularly in comparison to the
study of the direct problem. The earliest (engineering) papers on the subject were
published near the beginning of the nineteen-sixties. Mathematicians were thinking about
such problems much earlier; in 1923 Hadamard formalized the concept of well-posed, and

hence ill-posed, Lamm (1993). Due to the recent nature of the field, there are few

3
comprehensive books on the subject. Books by Kurpisz and Nowak (1995), Hensel

(1991), Beck et a1. (1985), and Beck and Arnold (1977) summarize and present various
methods and discuss practical applications of inverse methods. These aforementioned
books tend to be the broadest in coverage of inverse methods. Less comprehensive, but
more theoretical approaches are given in Alifanov et a1. (1996), Alifanov (1994), Ingham
and Yuan (1994), Murio (1993a), Baumeister (1987), Tikhonov and Arsenin (1977), and
Lawson (1974). These book present a more limited discussion, typically focusing on a
particular approach for inverse methods. They also tend to be more mathematically
rigorous. Additional references (most in Russian) are found in Kurpisz and Nowak (1995);
a relatively large Russian contingent is studying inverse problems.

As the importance and wide-spread application of inverse methods are realized, so too
have the demands on the complexity of the problems that can be solved. For twenty years
one-dimensional problems were the main subject of attention. Over the past decade multi-
dimensional problems have come under scrutiny, as well as inverse radiation problems and
most recently inverse convection problem. (Specific references are given in subsequent
chapters.) The demands on the complexity of the inverse problems that can be solved are
certainly increasing.

Multi-dimensional inverse problems are the main focus of this dissertation. In
particular, the estimation of coefficients in a describing partial differential equation and
the estimation of surface conditions from internal measurements for the two-dimensional
case are studied. These two problems were previously referred to as PE and IHCP,

respectively. The problems are closely related. In the former problem coefficients in the

4

describing partial differential equation are estimated. While the latter problem estimates
boundary conditions, or functions, in the describing partial differential equation.

An additional complexity addressed in this dissertation is the experimental applica-
tion of the multi-dimensional inverse methods. Although the contingency studying inverse
problems continues to grow, a relatively small faction applies the methods to experimental
studies, and even a smaller number have addressed multi-dimensional experimental appli-
cations. This dissertation studies inverse methods to experimentally estimate thermal
properties and surface conditions.

Experiments with carbon-carbon composite are used to study one and two-dimen-
sional inverse applications. Carbon-carbon composite, which is known to have thermo-
physical properties that depend on direction within the material (Loh and Beck, 1991), is
studied for the PE problem. For the IHCP a sequential method is proposed. This method is
hybrid in nature, although not completely, in that a sequential implementation on time of
the conjugate gradient method using an adjoint equation approach is proposed to solve the
IHCP. The sequential method was previously proposed for use with the function specifica-
tion method (Beck et al., 1985). Although a sequential implementation of the gradient
method has been proposed by Reinhardt and Hao (1996a), a very limited study of its
implementation was addressed for the one-dimensional problem. This is the first known
application of the method to a multi-dimensional problem. Issues of implementation and
benefits of such a sequential method are also investigated and discussed in detail.

Specifically there two major objectives, both related to multi-dimensional inverse
problems, of this dissertation:

1. Measurement of thermal properties of a carbon-carbon composite from
experiments with one- and two-dimensional heat flow. Two components of

5

thermal conductivity and the volumetric heat capacity are simultaneously
measured from transient temperature and heat flux measurements from

room temperature to 500° C.

2. Develop a conjugate gradient method using an adjoint equation approach
for solving the IHCP with a sequential-in-time concept. The sequential
concept is contrasted with the standard whole domain approach and varia-
tions of the sequential approach are proposed. Applications of these meth-
ods are studied with experimental data.

An outline of the dissertation is given. Carbon-carbon composite material and the
experimental aspects to estimate its thermophysical thermal properties are discussed in
Chapter 2 with the estimated one—dimensional properties presented. Chapter 3 gives the
two-dimensional thermal properties for the carbon-carbon and a comparison between the
one- and two-dimensional properties. A derivation of the conjugate gradient method with
adjoint equation for the IHCP is given in Chapter 4, including a sequential implementation
of the method. Chapters 5 and 6 investigate one- and two-dimensional applications of the

sequential method, respectively. A summary and conclusions are given in Chapter 7, and

recommended future work is suggested in Chapter 8.

Chapter 2

CARBON-CARBON THERMAL PROPERTIES: ONE-
DIMENSIONAL EXPERIMENTS

2-l.0 Introduction

Composite materials are defined as materials in which two or more constituent mate-
rials are combined to produce a resultant material that has different properties from the
individual constituents. Often the resulting properties (mechanical and/or thermal) of a
composite material are far better than that of the constituents, hence the advantage of the
composite materials. Carbon matrix carbon fiber composites are an important widely-used
family of such composite materials.

Carbon-carbon composites (CC) display the mechanical benefits of fiber-reinforced
materials, such as high strength-to-weight ratios, stiffness, and in-plane toughness, while
maintaining these mechanical properties at elevated temperatures. Carbon-carbon’s reten-
tion of mechanical properties at elevated temperatures (up to 2000 °C ) is unprecedented
(Savage, 1993). As a consequence, carbon-carbon is used in structural materials for space
vehicles, rocket nozzles, and aircraft brakes. Additionally, carbon-carbon is used for bio-
medical applications in hip replacement joints and for automotive applications in brakes,
clutches, engine blocks, and piston rings. Most automotive applications are associated
with motor sports, such as Formula I, or advanced applications. Presently, costs are too
prohibitive for commercial automotive applications.

6

7

I

The mechanical and thermal properties of the carbon-carbon composite are quite var-
ied. Depending on the weave and precursor used for the fiber and the method used to form
the matrix, the resulting properties (thermal and mechanical) will vary, possibly orders of
magnitude. In general, however, the properties tend to be anisotropic, varying with direc-
tion in the material. The complexity of the anisotropy depends on the weave and form of
the fiber used. For some carbon-carbon materials (such as the one studied herein) it is suf-
ficient to describe the properties based on two principal components, one describing the
properties in the direction parallel to the fibers and a second component for the direction
normal to the fibers. This type of material is referred to as orthotropic. Fortunately in the
material studied, the fibers are oriented to align the fiber direction to coincide with the
planes of the outer surfaces.

Although carbon-carbon can maintain structural integrity at extreme temperatures, it
does so only in an inert atmosphere. At moderate temperatures (500 °C ) in an air atmo-
sphere oxidation will occur (Bines, 1993). To prevent oxidation and extend the tempera-
ture range the carbon-carbon requires protection. Typically, oxidation protection is
obtained by coating the surface with a form of silicon (silicon-carbide, SiC) through one
of many application processes (chemical vapor deposition (CVD), infiltration (CVI), or
reaction (CVR), or “pack coated”). After coating the surface, a film (tetraethylorthosili-
cate, TEOS) is applied to the oxidation protective coating to seal the coating. The TEOS
film significantly decreases relative weight loss (oxidation) through the porous SiC coat-
ing. Since the SiC has a relatively high thermal conductivity (Incropera and Dewitt, 1990)

and appreciable thickness, the properties of the carbon-carbon are even more complex

8

with the oxidation protective coating. (For this investigation, the effects of the CC and SiC
are lumped and effective thermal properties are estimated.)

The carbon-carbon composite is known to be a thermally complex material. The sim-
plest description for the thermal properties of the material, beyond isotropic, are orthotro-
pic. Hence, multi-dimensional experiments are required to fully characterize the
properties. (A series of one-dimensional experiments would not verify the adequacy of the
orthotropic assumption.) Furthermore, with the number of different fiber/matrix variations
and manufacturing methods a wide range of thermal properties are possible (Savage,
1993). Consequently, methods, experimental and analysis, are required to measure the
multi-dimensional thermal properties for accurate modeling of the material. Parameter
estimations techniques permit the measurement of thermal properties from appropriate
experimental measurements. These methods are applied and studied to estimate the ortho-
tropic thermal conductivity and isotropic volumetric heat capacity of the carbon-carbon.

The remainder of this section gives a problem description and literature review.
Experimental aspects are discussed in Section 2-2.0. Parameter estimation methods for
analyzing the experiments are covered in Section 2-3.0. Section 2-4.0 presents and dis-
cusses the thermal properties estimated for the one-dimensional experiments. Experimen-

tal uncertainty is analyzed in Section 2—5.0.

2-1.l Problem Description

The development of composites materials is proceeding at an ever increasing rate.
Advanced materials are used in automotive, airplane, and aerospace applications. The
strength-to-weight ratio and survivability in harsh conditions are two main advantages of

these materials. However, the emergence of these materials requires experimental

9

techniques and solution methods to determine their thermal properties. This chapter is
about the estimation of thermal properties from temperature and heat flux measurements
for a carbon-carbon composite material.

The investigated material is an advanced carbon matrix-carbon fiber material made by
Carbon-Carbon Advanced Technologies, Inc. of Fort Worth, Texas. There are six speci-
mens, 7.62cm square and 0.953cm thick. (Testing, presented herein, was performed on
two of these specimens.) The specimens are described (by the manufacturer) as CCl 2-D
composite made from fiberite K-641 fully densified, SiC pack coated with sealant.
Because the specimens were not flat as received, the specimens were ground to produce

flat surfaces.

2-1.2 Literature review

Several methods have been used to measure the thermal properties of composite
materials. These methods include: the guarded hot plate (Lee and Taylor, 1975), flash
method (Lee and Taylor, 1975, Taylor et al., 1985, and, Harris et al., 1982, and parameter
estimation (Beck and Osman, 1991 and Gamier et al., 1992). In addition, modeling of the
composite has been investigated to derive properties from constituent or laminate proper-
ties; see Kulkami and Brady (1997), Balageas and Luc (1986), Han and Cosner (1981),
and Chamis (1973). In comparison to the other methods, parameter estimation provides
interaction between the model (analysis) and the experiment. For example, in the flash
method a small sample is heated on one side for a very short duration, typically with a
laser, and the temperature is measured on the opposite side. Early methods to calculate the

thermal diffusivity (0:) used the relationship

10

2
a=i=91fli (2-1)

PC ’ 1/2
where (1,2 is the time required for the temperature to reach one-half the maximum tem-

perature rise and l is the sample thickness. An advantage of this method is that it is rela—

tively simple. However, the method has certain disadvantages:

The (heat conduction) model equation is not necessarily satisfied nor checked

Typically few measurements are used

The properties k and pC cannot be obtained independently

PP’P?‘

Generalizing these methods, i.e. multi—dimensional, more complex model, or temper-
ature dependent thermal properties, is difficult, if not impossible

5. Use of statistics to quantify estimates or improve the experiment is not easily done

Parameter estimation does not posses the limitations of this other method. Beck (1996)
discusses applying parameter estimation techniques to the flash diffusivity experiment.
Beck shows by observing parameter estimation results that the model can be refined to
improve the accuracy of the estimated thermal diffusivity. Not only by improving the anal-
ysis, but by using the analysis to better understand the experiment. Using equation (2—1) to
calculate the thermal diffusivity does not permit such an analysis or provide insight.
Hence, the more general method applying parameter estimation techniques is used in this
dissertation.

The theory of the techniques in this research using electric heaters for estimating the
thermal properties is detailed in the book Beck and Arnold (1977, see Chapter 7 in partic-
ular). A consequence of using electric heaters is that typically both thermal properties can

be estimated simultaneously. This research uses the electric heaters.

11

Other methods with the known heat flux and transient measurements are described in
Beck, et al. (1991), Scott and Beck (1992a,b), and Beck and Osman (1991). The first of
these papers uses an internal heat flux transducer. The papers by Scott and Beck (1992a,
1992b) relate to composite materials during and after cure. A method to sequentially esti-
mate thermal properties by mathematically connecting a series of discrete experiments
with various temperature ranges is presented by Beck and Osman (1991). Loh (1991) has
given results for experimentally determining two components of the thermal conductivity
for a carbon-carbon composite. The papers by Gamier et a1. (1992,1993) describe a
method for estimating thermal properties without requiring the temperature sensor inside
the specimen(s). This was also attempted in the present work, but the method (at least for
one-dimensional cases) is only appropriate for materials with relatively low thermal con-
ductivities because of the contact conductance. The carbon-carbon composite material
does not have low values, however. Also relevant is the study of optimal experiments,

which are discussed in Chapter 8 of Beck and Arnold (1977) and in Taktak et a1. (1993).

2-2.0 Experimental Aspects

2-2.1 Experimental Description

A sketch of the experimental set-up used to estimate the thermal properties of the car-
bon-carbon material is shown in Figure 2-1. It consist of two nominally identical carbon-
carbon specimens (7.62cm x 7.62cm x 0.914cm) and ceramic insulations (7.62cm x
7.62cm x 2.54cm, Zircar Products Inc., Florida, NY) with a mica heater assembly (Ther-

mal Circuits, Inc., Salem, MA., (2(T = 33 Ohms ) located between the identical

room)

halves. Five thermocouples (Type E, 0.254mm nominal wire diameter) are embedded on

12

ceramic insulation

 

084

 

 

 

 

 

 

('034) 13 composite specimen 14
j 8 9 10 11 12
'l/IEIIIIIIILn/IIIAXXXELXxxxxl—I.xxxxxxxxxxxxxElxxxx
T 3 4 5 L5] .91
/ 1 composite specimen 2 V (.36)
mica
heatelr;l
assem y . . . 2.54
ceranuc insulation (1.00)
1L.
L 7.62 (3.00) _

 

_ '1

All dimensions in cm (in.)
(nor DRAW TO SCALE)
0 thermocouple Ad"? heat"
82 Inactive heater (2D)

Figure 2-1 Schematic of experimental set—up used for estimating thermal properties of
carbon-carbon composite. For one-dimensional experiments all heaters are activated.

13

the surface of each carbon-carbon specimen at the heater/specimen interface. The thermo-
couples (insulation removed) are attached with electrically insulating high temperature
cement into grooves (0.38 mm by 0.46 mm) that extend the length of the specimen. Two
thermocouples are located at each interface of the carbon-carbon specimen and the
ceramic insulation. The entire set-up is mounted between two 3.18 mm thick aluminum
plates that are connected with threaded rods and hold the numerous layers firmly in place;
the apparatus is placed in a furnace, which allows variation of the initial temperature. Fur-
ther details of the experimental procedures are discussed in Ulbrich et a1. (1993).

The experiments are conducted and processed using a 12-bit data acquisition system
(National Instruments) with a 486 PC. The system provides accurate data acquisition with
minimum sampling intervals in the usec range. Two eight-channel data acquisition boards
are linked providing sixteen channels of data acquisition. The system controls and
acquires the power (voltage and current) delivered to the heaters and acquires the therrno-
couple voltages. The current is acquired by measuring the voltage across a known resis-
tance. The heat flux is calculated from the power measurements assuming the heating is
uniform over the heating area (7.62 cm x 7.62 cm) and divides equally to the symmetric
experimental halves.

The measured temperatures are averaged on opposite sides of the heater assembly to
determine the temperature at each location. The location of the thermocouples are shown
in Figure 2-1 and given in Table 2-1. The sensors that are embedded in the specimen are
assumed to measure the temperature at the surface of the specimen. Since thermal conduc-
tivity normal to the fiber is much greater than the conductivity of the insulation

(k; cc » kins) , small temperature gradients near the specimen / insulation interface. The

14

Table 2-1 Sensor locations in experimental apparatus

 

Location cm (in)

 

Sensor Sensor

Location Number x y
A 3,8 0.89 (.35) 0
B 4,9 1.91 (.75) 0
C 5,10 3.18 (1.25) 0
D 6,11 4.45 (1.75) O
E 7,12 6.73 (2.65) 0
F 1,13 1.27 (0.5) 0.91 (0.36) (Ly)
G 2,14 6.35 (2.5) 0.91 (0.36) (Ly)

non-embedded sensors are assumed to measure the temperature at the backside of the car-

bon-carbon specimen.

The thermal experimental model for the one-dimensional case is shown in Figure
2-2. All outer surfaces are assumed to be insulated, except for the surface where the
energy is introduced by the heater. The energy to the heater is assumed to divide equally
between the two halves and emanate from the middle of the heater assembly
(y = —0.042 cm). The adequacy of the assumed insulated outer surfaces for the model can
be verified by a comparison of heat losses from natural convection with the anticipated
applied heat flux from the heater assembly. Since a temperature rise of 20 to 25°C above
the ambient is expected for a typical experiment, the heat loss is mainly due to natural con-
vection (h z 4W/m7-C). These losses are negligible (qwnv z 100W/m2) in comparison
to the applied heat flux, which was a minimum of 7800W/m2 for the one-dimensional

testing. Hence, insulated boundary conditions on the outer surfaces are applicable.

 

 

q = 0
11"
1’2 mica ceramic insulation L-
heater ms
assembly
11
. 1,2,13,14. 1
y composrte specrmen Ly
3,4,5,6,7,8,9,10,11,12
'l [A] III ll’ll'l

 

f/I’I I 1

11/111111it1111111111111111111111117111'1111'1111' 1111111111111 f
L

 

q(t)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mica

Figure 2-2 Heat transfer model for one-dimensional experiment

16

The thermal model is mathematically represented as three coupled problems. In
the analysis of a particular experiment the thermal properties are assumed constant in the
model. Temperature dependence of the thermal properties is developed by combining
experiments over the temperature range. The first problem is for the mica heater/contact

conductance

2
Tmi 8T ic
Strung; C0 = (prnicagt'm a’_Lmica<y<0’t>O (2'2)
8T .
42...?!“ = qm (2-3.,
y y = 0
Tmica(y,0) = To, —LmicaSySO (2-3b)
The second problem is for the carbon-carbon composite
ke 66— Ce 8T“ 0 L t>O 24
y,ch _(p )CC-g; ’ <y< y’ (')
ch(y,0) = To , OSySL), (2-5)
The third problem represents the insulation material
e ins e aTins
kins-3:23 = (901,255; .Ly<y<Ly+Lins,t>0 (2-6)
aTins
_a_, = 0 (2-7a)
} r=Lm+A
Tim.(y, 0) = To, L), S y _<. Ly + Li,” (2-7b)

The problems are coupled through interface conditions that assume perfect contact

ke arm,m(o,z)_ e BTCC(O,t)

mica av _ .y, CC ay 9 (2'83)

17

Tmica(0’ t) = ch(0’ t) (2-8b)

6T," (L ,t) aT (L ,t)
kins Say y = k; cc CCay)’ (2-8C)
Tins(Ly’ t) = ch(Ly’ I) (2-8d)

Although perfect contact is not reasonable between the mica heater and carbon-carbon,
effective properties are used to represent the mica heater and contact conductance due to
imperfect contact. Effective thermal conductivity is defined

8T

ke : —/— -
q a)? (2 9)

_ . T . .
where q rs the average heat flux and 9— 13 the average temperature gradient. Both aver-

3y

ages are in the direction normal to the heat flow. The effective volumetric heat capacity is

defined

(pC)"’ = ‘l/IpCdv (2-10)
v

where V is the volume. By experimentally measuring effective properties for the mica
heater assembly and insulation the contact conductance due to imperfect contact is
accounted for in the model. In addition, the effective properties of the mica heater account
for the effects of other materials in the heater and the cement used to install the thermo-
couples in grooves in the carbon-carbon composite. Modeling the effects of several mate-
rials with an effective property is much easier (and in most cases more accurate) then
accounting for all effects individually. Effective properties of the carbon-carbon represent
the non-homogenous construction of the material. To protect from oxidation a thin layer

(approximately 10% of the overall thickness) of silicon-carbide protects the outer surfaces

18

of carbon-carbon. Consequently, effective properties that represent the carbon-carbon and

silicon-carbide are estimated.

Including the mica heater and ceramic insulation in the model (Figure 2-2)
requires the thermal properties of the materials to be known, or also estimated, to deter-
mine the properties of the carbon-carbon specimen. Neglecting the mica heater in the
model is not appropriate because the contact conductance results in a large temperature
drop across the heater. If the heater is neglected, the carbon-carbon properties will incor-
rectly reflect this effect. Also, including the heat loss to the insulation, instead of assuming
a perfect insulated condition (at the specimen/insulation interface), increases the accuracy
of the properties estimated for the carbon-carbon. If known (tabulated or published) prop-
erties are used for these materials (mica and insulation), several problems arise. First, ther-
mal properties are typically not known very accurately. Second, contact conductance
between adjacent layers is typically not negligible and must be considered. Third, the
ceramic insulation was treated by spraying with a ridigizing material, possibly changing
its thermal properties. By experimentally estimating effective properties of these materi-
als, these problems are less influential in the estimated properties of the carbon-carbon.

This approach requires additional experimental work, however. The experimental
conditions (length of experiment and heating duration) to determine the properties of the
mica and insulation are quite different from the conditions necessary to estimate the prop-
erties of the carbon-carbon. A separate series of tests are performed, one set with the car-
bon-carbon specimen removed, to determine the properties of the insulation and mica
heater assembly. These are effective properties, which will account for any contact con-

ductance due to imperfect contact between the different material layers.

19

With the separate series of experiments the effective thermal properties of the ceramic
insulation and mica heater assembly were determined. However, when the set-up is recon-
figured to test the carbon-carbon composite, the contact conductance may vary. Therefore,
for each experiment to determine the properties of the carbon-carbon composite, a short
duration experiment is conducted to characterize the effective thermal properties of the
mica heater and contact conductance. The short duration experiment is conducted such
that the thermal properties of the mica heater, including the contact resistance, are impor-
tant in the thermal model, but the properties of the carbon-carbon are not important.
Importance of the materials in the model is quantified by the sensitivity coefficients (dis-
cussed in Section 23.0). A duration is selected that is sensitive to the properties of the
mica, but relatively insensitive to the properties of the properties of the carbon-carbon.

For a typical set of experiments to determine the properties of the carbon-carbon
composite, the furnace is set at a given temperature and the experimental apparatus is
allowed to reach a uniform temperature (usually this required several hours). A short dura-
tion experiment (heating for approximately 2 seconds) is run first, which is used to esti-
mate the effective properties of the mica heater. After allowing the set-up to stabilize at a
uniform temperature, a second experiment is run. The second experiment (heating for
approximately 20 seconds) is longer in duration and is used to estimate the properties of

the carbon-carbon composite specimen.

2-2.2 Experimental Design
During the initial stages of this research a different heater design was used. The heater
was constructed of a Kapton material with resistance temperature sensors (RTD) integral

in the heater assembly. This approach had the advantage of being non-intrusive and

20

required no machining to install thermocouples. A major difficultly with this design was
that the magnitude of the contact conductance resulted in the one-dimensional thermal
resistance of the heater/contact conductance being on the same order of magnitude as the

composite specimen in the thermal model

L L
y z kapton (2_11)

e e
k J” CC kkapton

where kzapton is the effective thermal conductivity of the kapton heater/contact conduc-
tance and Lkapmn is its thickness. Hence, the measured RTD temperature depends as
much on the thermal properties specified for the heater/contact conductance as it does on
the thermal properties specified for the composite specimen. (Sensitivity coefficients, dis-
cussed in Section 2-3.0, for the two materials would be comparable in magnitude.) In gen-
era], when other materials are present in the model, their effect on the temperature should
be small in comparison to the material for which the properties are sought. For this case,

the sensitivity coefficients (discussed below) for k: C should be larger, compared to the

c
sensitivity coefficients for kzapwn . To improve the accuracy and reduce the dependence on
the thermal properties of the heater/contact conductance, the sensors were embedded in
the surface of the specimen.

The choice of the boundary condition on the backside of the specimen is dependent
on the thermal conductivity of the test specimen. For a test specimen with a relatively high
thermal conductivity (which this carbon-carbon does have) it is more practical to simulate
an insulation condition. An insulating boundary condition is less sensitive to the contact

conductance at this surface. However, for a relatively low thermal conductivity material it

is more practical to simulate a temperature boundary condition. For a test specimen with a

21

relatively low thermal conductivity, the boundary temperature is again less sensitive to the
contact conductance. Based on optimality criteria for the one-dimensional case with an
applied surface heat flux, a specified temperature boundary condition on the back surface
is an optimal experiment for estimating thermal conductivity. An insulating condition is
optimal for estimating the volumetric heat capacity, Beck and Arnold (1977).

An approximate insulation boundary condition on the back of the specimen is chosen,
instead of a temperature boundary condition (which was initially used), to minimize the
sensitivity to the contact conductance at this location. With the insulation boundary condi-
tion, the heat flux (and temperature gradients) at the specimen/insulation interface is small
and the temperature is not sensitive to the contact conductance or location for the non-
embedded sensors. A temperature boundary condition, however, has a larger heat flux at
the boundary and the temperature (at the sensor) is very sensitive to the magnitude of the

contact conductance and the specified location for the sensors at the interface.

2-3.0 Analysis Procedures

The techniques to estimate thermal properties are detailed in a book by Beck and

Arnold (1977). The basic process involves minimizing a sum-of-squares function

5 = (Y—T)TW(Y—T)+(u—b)TU(u—b) (2-12)

where Y and T are vectors of the measured and calculated temperatures and W is a
weighting matrix (typically the identity matrix). The last term in equation (2-12) serves as
regularization or allows for'the inclusion of prior information about the thermal properties.
It contains the difference between the prior information u and present estimates b with a

symmetric weighting matrix U. To determine the thermal properties the function S is

22

minimized with respect to the thermal properties, i.e., 171 = (pC):c, b2 = k; cc, and

b3 = ke This is accomplished by setting the first derivative of S with respect to each

x,cc'

parameter equal to zero, and solving for the estimated parameters (3). The resulting

expression for the estimated thermal properties (Beck and Arnold, equation 7.4.6, 1977) is

x(k+l) _

A k A A
b b( )+ P(k)[XTW(Y — T) + U(u — b(k))] (2-13)

P“) = [XTWX+ U1“ (2-14)

The superscript (k) defines the iteration number, iteration is required even for the linear
conduction problem, due to the non-linear nature of the estimation problem; the sensitivity
coefficients X depend on the parameters (thermal properties). The columns of the sensi-
tivity matrix are the first derivative of temperature (for each time and sensor location) with

respect to the parameters

X = [X19], sz, . . . ,Xbp] (2'15)
.. 61
" ab.

1

X b (2-16)
The sensitivity coefficients can provide considerable insight to the estimation problem and
aid in the design of the experiment for optimum accuracy in the estimates (Beck and

Arnold, Chapter 8, 1977). One criteria for an “optimal” experiment, valid for additive,

zero mean normal errors in Y, and errorless independent variables, is to maximize

A = IXTXI (2-17)

This criteria is appropriate because it corresponds to minimizing the volume of the
confidence region for the estimated parameters. By studying the experimental design prior
to conducting experiments, such that equation (2-17) is maximized, the maximum

information is available from an experiment. In addition to the criteria proposed in

23

equation (2-17), a constraint requiring a fixed number of observations and/or the same
temperature rise, may be required to provide consistency in comparing experimental
designs with multiple sensors. Taktak et al., (1993) discuss the design of an experiment to
estimate thermal properties of composites with relatively low thermal conductivity.
Additional investigations on experimental design concerning the optimal placement of
sensors (Fadale et al., 1995a) and the design of experiments using uncertainty information
(Fadale et al., 1995b and Emery and Fadale, 1996) use the Fisher information matrix.

The sensitivity matrix and the temperatures are linearized about the thermal proper-
ties from the previous iteration in equation (2-13) and (2-14). Iteration continues until con-

kl Ak ak
(+)_b() 3 8b“,

vergence of the estimated parameters is reached, as defined by, |5
where e is a small number to quantify convergence, such as 8 = 0.0001 .

These solution procedures are implemented with the computer program PROPlD
to estimate the thermal properties. PROPID provides a means to estimate thermal proper-
ties of multi-layer bodies from appropriate measurements. Thermal conductivity and volu-
metric heat capacity may be determined simultaneously and for more than one material, if
desired. Layers of different materials may also be lumped and effective thermal properties

determined. Gamier et a1. (1992) performed experiments on materials with well-known

published thermal properties to support the accuracy of PROPID.

2-4.0 Results and Discussion

Three types of experiments were conducted. The first type removed the carbon-car-
bon from apparatus and replaced it with a 1.25 cm thick insulation. These type-one exper-

iments estimated effective thermal properties of the mica heater and insulation.

24

Experiments of type-two and type-three replaced the carbon-carbon specimen in the
experiment. Type-two experiments heated for a short duration, providing information
about the effective properties of the mica heater. Type-three experiments estimate the
properties of the carbon-carbon. Results from experiments of type-one and type-two are
discussed in Section 2-4.1 and type-three is discussed in Section 2-4.2.

The reason that three types of experiments were used, is to quantify three separate,
quite different, effects. Type-one experiments, with the carbon-carbon composite
removed, provide effective thermal properties of the insulation. Experiment type-two was
relatively short in duration, so that temperature is mainly influenced by the properties of
mica heater, contact conductance, and associate materials used to install thermocouples.
Type-three experiments were longer in duration and characterized the properties of the
carbon-carbon composite. Each of the three experiments focuses on an individual aspect

of the thermal model and has minimal effects from other material’s thermal properties.

2-4.1 Effective Properties for Mica Heater and Insulation

An independent set of experiments was conducted to determine the effective thermal
properties of the insulation material (Ulbrich, 1993). These tests used a set-up similar to
Figure 2-1, except that the carbon-carbon specimen was replaced by a 1.25 cm-thick piece
of insulation. The goal of these experiments was to estimate the thermal properties of the
insulation. The results are given in Table 2-2. The values estimated for the insulation are
considerably higher than the values reported by the manufacturer. This is not a surprising
outcome. Manufacturers do not typically measure thermal properties. Also, this material
was treated with a solution to make it more structurally strong. The effective properties for

the mica heater were also measured. Because these values are sensitivity to the contact

25

Table 2-2 Thermal properties of the ceramic insulation

 

 

 

 

 

 

 

 

 

 

 

Present Investigation Manufacturer
Temperature kins (p C )fnsx 10-6 kins
°C W/(m°C) J/(m3°C) W/(m°C)
40 .088 .419
175 .093 .570
200 .055
600 .1 10

 

 

Table 2-3 Effective thermal properties for the mica heater assembly

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Initial . kin-ca (pC)fm.wx10"6
Experimental Temp 0 3
Case (°C) (°C) W/(m°C) J/(m oC)
1010@3(),1 30 0.066 0.142 +/- 0.002 2.030
1010@41_1 41 0.070 0.131 +/- 0.002 2.030
1007@134,1 134 0.067 0.102 +/- 0.001 2.030
1007@143,1 143 0.072 0.1 10 +/- 0.002 2.030
1008@180,1 180 0.057 0.125 +/- 0.001 2.030
1012@195,1 195 0.053 0.123 +/- 0.001 2.030
1008@245.1 245 0.070 0.1 14 +/- 0.001 2.030
1013@254,1 254 0.057 0.123 +/- 0.002 2.030
1008@295,1 295 0.058 0.123 +/- 0.002 2.030

 

 

26

conductance they were measured again after the apparatus was re-assembled with the car-
bon-carbon specimen. These experiments were referred to as type-one experiments.

Having established the properties of the insulation, type-two experiments were con-
ducted to determine the effective properties of the mica heater. It was found that estimat-
ing the thermal conductivity and volumetric heat capacity simultaneously is not possible
due to correlated sensitivity coefficients (see Beck and Arnold, 1977). Therefore, the volu-
metric heat capacity of the mica heater is set at the value estimated during the experiments
to determine the properties of the ceramic insulation (type-one). The effective thermal
conductivity of the mica heater is estimated. The thermal conductivity is estimated
because it is most influenced by the contact conductance, which may have changed from
the previous estimates.

The effective thermal conductivity estimated for the mica heater assembly, assuming
(pC);,-ca = 2.0><106 J/m3C, demonstrated no specific trend with temperature over the

range (30—295)°C . See Table 2-3. The largest value is ke (30°C) 2 0.14 W/mC

mica

e
mica

and the smallest is k (134°) = 0.10 W/mC . For temperatures up to 295°C the ther-
mal conductivity of the mica heater is estimated at each temperature. After which, the
magnitudes of these thermal properties were shown to minimally affect the outcome of the
estimation for the carbon-carbon properties. A variation of 50% in the thermal properties
specified for the mica heater results in variations in the estimated properties of the carbon-
carbon that are within the magnitude of the confidence intervals (discussed below). For
temperatures greater than 295°C , the thermal conductivity estimate of the mica heater at

295°C is used in the analysis. The fact that it is difficutlt to estimate the thermal proper-

ties of the mica heater suggests that these properties are not significant in the estimation of

27

Table 2-4 Thermal properties estimated for the carbon-carbon composite from
one-dimensional experiments and analysis

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. . e
Experimental 113:2: 6 ky’ CC ( p C):c3x 10-6
Case (°C) (°C) W/(m°C) J/(m °C)

1010#30.1 31 0.128 3.40 +/- 0.05 1.42 +/- 0.01
1010#41.1 42 0.114 3.48 +/- 0.04 1.47 +/- 0.01
1012#94.1 95 0.161 3.85 +/— 0.07 1.74 +/- 0.02
1012#108.l 109 0.161 3.93 +/- 0.07 1.81 +/- 0.02
1011#143.l 143 0.121 4.12 +/- 0.07 1.90 +/- 0.02
1011#159.1 159 0.089 4.22 +/- 0.06 2.08 +/- 0.01
1012#195.l 195 0.091 4.35 +/- 0.06 2.20 +/- 0.02
1013#259.1 259 0.075 4.60 +/- 0.04 2.47 +/- 0.02
1008#295.l 295 0.090 4.76 +/- 0.05 2.52 +/- 0.02
1013#304.1 304 0.103 4.74 +/— 0.05 2.58 +/- 0.02
1023#403.2 403 0.100 4.86 +/- 0.05 2.76 +/- 0.02
1023#455.2 455 0.082 4.93 +/- 0.03 2.88 +/- 0.02
1023#508.2 508 0.096 4.99 +/- 0.05 2.97 +/- 0.02
1023#571.2 571 0.278 3.93 +/- 0.23 3.06 +/— 0.10
1023#623.2 623 0.205 3.73 +/- 0.15 3.23 +/- 0.07

 

 

the carbon-carbon properties (or thermal model). This is a desirable characteristic for
additional materials (i.e. mica heater and ceramic insulation) in the experimental model. A

similar insensitivity is demonstrated for the thermal properties of the insulation.

2-4.2 Effective Properties of Carbon-Carbon
The properties estimated for the carbon-carbon composite are given in Table 2-4. The

first column identifies the experimental case and the second column the initial

28

temperature. The third column is the estimated standard deviation of the analysis, which is

defined as

1 I .I 2 0.5
O = [W2 2 (Yij_ 7.0)] (2-18)

1 = l j = l
where N t is the total number measurements (from both sensor locations, N t = I + J) and

N p is the number of estimated parameters. Variables I and J are the number of measure-
ment times and number of sensors, respectively. The last two columns present the esti-
mated properties with the associated confidence interval (calculated by PROPID and
discussed below).

Insight about the estimated properties is gained by plotting the properties as a func-
tion of temperature; Figure 2-3 gives a plot of the one-dimensional properties of the car-
bon-carbon composite as a function of the temperature, with the estimated confidence
intervals. Note that the analysis assumes that the thermal properties were constant during
an experiment, but varied between experiments; the initial temperature of the experiment
is used to plot the properties. Using an F-test (Beck and Arnold, chapter 6, 1977) it is con-

cluded that a second order (in temperature) model adequately represents the results. The

equations determined for the properties with a least squares fit are

 

 

 

 

k" 342 338 T40 186 T-TO 2 9
= . . — . 21
y, CC + (T3 ‘ T0) (T3 _ To) ( )
C" 141 265 T40 111T—TO 2 106 0
(P )CC — [. + . (T3—TO)_ . (T3—TOHX (22)

where T0 = 31°C and T3 = 508°C are the minimum and maximum temperatures. By
normalizing the temperature in equation (2-19) and (2-20) the coefficients in the

describing equations have units equivalent to the respective thermal properties. A physical

appreciation of the magnitudes is more easily realized in this form. The relationships are
also shown in Figure 2-3. The thermal conductivities determined from the last two
experiments are not used in the least squares analysis. The validity of the thermal
properties estimated at temperatures of 571°C and 623°C is uncertain. The residuals and
confidence intervals for these two cases are relatively large. When the set-up was
disassembled after these experiments, the outer coating of the specimen (which protects
from oxidation) had separated from the inner composite material. Because the estimated

thermal conductivities are lower it is quite possible that the failure of the specimen

29

occurred just before or during these experiments.

Thermal Conductivity [W/(mC)]

 

 

 

 

0 '100'200'300'400'500'600

Temperature (C)

10.01 . . . r . 3.5E+06
9.0- 5

- c " —3.0E+06
8.0- (p )“ .
7.0.. -2.5E+06
50': kg,“ -2.0E+06
5.0- .
450‘ :1.5E+06
30, -1.0E+06
2'0- .501: 05
1.0- _ ' +
0.0 0.0E+00

Figure 2-3 Temperature dependence of thermal properties estimated from one-

dimensional experiments

Volumetric Heat Capacity [J/(m3C)]

30

The details of the parameter estimation are quite similar at different initial tempera-
tures. For this reason, and for brevity, only one case (1023#508.2), at an initial tempera-
ture of 508°C , is discussed. The experimental data for this case are presented in Figure 2-
4. The temperatures shown represent the average of all thermocouples at the referenced
locations. The results closely approximate the case of a finite slab heated with a constant
heat flux at the surface and insulated at the backside. Initially, the temperature at the
heated surface increases rapidly, while the backside temperature remains constant. Later
both the temperatures at the heated surface and backside increase linearly with time until
the power to the heater is turned off. Then the temperatures at both surfaces tend to
approach the same temperature, demonstrating that there is little heat loss to the insulation

(even though it is considered in the analysis).

 

 

 

 

 

 

525.0
520.0-
8 515.0-
a
a
e
g» 5100- .42
o E
i” E
505.0- Heat Flux {15000.0 3
_-1oooo.0 a
35000.0 ‘5};
500.0 - . . . 4 4 . , 0.0 m
0 10 20 30

time (seconds)

Figure 24 Experimental data for one-dimensional case 1023#508.2

31

In addition to estimating the thermal properties, PROPlD provides some means to
quantify the accuracy of the estimates. The previously discussed estimated standard devia-
tion, equation (2-15), provides an indication of how well the calculated match the mea-
sured temperatures. The magnitude of ("I can be compared to the temperature rise of the
experiment, which is approximately 15°C . Except for the last two experiments, 6 is
within 1% of the temperature rise. Also given with the estimates of the parameters is a
confidence interval (Beck and Arnold, Chapter 6, 1977). The calculation of the confidence
intervals has associated assumptions. First, the model for the experiment is correct. Sec-
ond, the dominant errors in the analysis are in the temperature measurements, modeled
with a first order auto-regressive model, and the errors are not biased (Beck and Arnold,
1977).

Other quantities can be observed to demonstrate the accuracy of the estimated proper-
ties. These include the sequential estimates of the properties, the residuals, and the sensi-
tivity coefficients. The quantities are important to provide insight to the estimation, as well
as insight to the experiment. Observing them can help improve the experiment and support
the accuracy of the estimated properties. Each is discussed below.

The sequential estimates demonstrate how the estimated properties vary as additional
measurements are considered. The analysis assumes that the estimates are linearized about
the converged parameter values using all experimental data. Figure 2-5 shows the sequen-
tial estimates for this case. The sequentially estimated property, at time I, , represents the
outcome if only data up to that time is used in the analysis (and linearized about the con-
verged property values). In other words, if the data is analyzed by adding one data pair,

Y( y = 0) and Y( y = Ly) , at each time, it shows how the estimated properties change as

32

 

 

 

 

6.0 . . . . . . . e 3.2E+06
:Q —\V’i (pC)CC :
Lg A e -2.8E+O6 (,9
i 5.0" V— ky, cc g
E L2.4E+00 a
.§ 4.0- r?
,2 -2.0E+06 g
g 8*
g 3.0- -1.65+06 3
O a
U 0
"a -1.2E+06 f
E 2.0- '5
2 -a.01:+05 “é
{-1 :r
1.0- .5
-4.0£+05 >

0.0 . r 0.OE+00

 

 

0 5 1'0 1'5 2'0 2'5 30 :55 40
time (seconds)

Figure 2-5 Sequential parameter estimates for one-dimensional case 1023#5082

an additional data pair is added to the analysis. Initially the sequential estimates vary
because there is not enough information to determine the parameters. However, as more
data is considered, the property estimates approach constants and the linearization approx—
imation is accurate. Meaning, if the experiment (or analysis) is ended at 15 seconds, the
estimated properties would not differ significantly from the properties at 30 seconds. In
general, for a good estimation, the sequential estimates converge to a constant and are
fairly steady with time. For times greater than fifteen seconds, the sequential estimates of
k; CC and (pC):c vary 0.9 and 3.5%, respectively. These values can be compared with the
confidence regions predicted in Table 2-4 for this case of 1% and 0.67%, respectively.
These sequential results for (pC):c at 508°C are not as accurate as the confidence inter-

val in Table 2.4. Uncertainty in the other experimental measurements account for the dif-

ference; the confidence intervals include error in the temperature measurements only.

33

The residuals are related to the estimated standard deviation and are calculated using

and represent the difference between the measured and calculated temperature for a partic-
ular time (ti) and sensor location (j) . The estimated standard deviation gives an indica-
tion of the magnitude of the residuals; the signs and magnitudes of the residuals can
provide considerable insight. Figure 2-6 presents the residuals for this case. The magni-
tude of the residuals is approximately 0.1°C . The residuals are correlated; most of the
residuals are positive during the heating. This outcome may signify that some inconsis-
tency exists in the model or that a small effect was omitted. However, the magnitudes of
these residuals are small, within 1% of the temperature rise during the experiment, indicat-

ing that errors in the model are minimal.

 

 

 

 

0.5 l I I 1 l 1 I
0.3- —
. .4
- '01 . -"' -
A :I. ‘ l." V‘ \ ' .
U 0°1- . I I "5 '1“: 1i -
v II. It. .5 I:
a .l I ~ I . WI. 1‘ q
E t . 2'1: 3'” t r
U) _ _ U s. c 3 . ' _
32 0.1 I; 615...": .1
-0.3- —— sensor (y = 0) _
----- sensor (y = Ly)
‘005 l l l r I I

051015 20 25 30 3540

time (seconds)

Figure 2-6 Temperature residuals for one-dimensional case 1023#508.2

34

Sensitivity coefficients are the first derivative of the temperature with respect to the
parameters, thermal conductivity and volumetric heat capacity. They are indicators of how
well designed the experiment is. In general, the sensitivity coefficients are desired to be
large and uncorrelated (linearly independent). A sense of the magnitude of the sensitivity
coefficients is gained through normalizing the sensitivity coefficients by multiplying by
parameters, resulting in units of temperature for the normalized sensitivity coefficients. A
comparison is then possible with the temperature rise of the experiment. For a well
designed experiment, with boundary conditions similar to the case investigated in this
study, the sum of the normalized sensitivity coefficients for the thermal conductivity and
volumetric heat capacity is nearly equal to the negative of the temperature rise (equal only
if perfectly insulated at y = Ly ). Sensitivity coefficients are useful in the design of exper-
iments, i.e., determining the heating and experiment duration, location of sensors, heating
area (for two-dimensional case), etc. A study of the sensitivity coefficients, prior to per-
forming experiments, can lead to better experiment designs.

Figure 2-7 shows the sensitivity coefficients for the representative experimental case
(1023#508.2). The sensitivity to the thermal conductivity and volumetric heat capacity are
shown for both sensor locations. The magnitudes are about equal to the temperature rises,
which is a good feature. Notice that the sensitivity coefficients are correlated (linearly
dependent) for times up to 10 seconds for the sensor at the surface of the specimen
(y = 0). This is similar to the situation that resulted in only being able to estimate one
parameter for the effective properties of the mica heater in the analysis of the short dura-
tion experiment. In this case, however, information is available from another sensor

. . . . . e e -
(y = Ly) where the normallzed sens1t1vrty coefficrents for k“ C and (pC)CC have qurte

C

35

different shapes (i.e., not correlated). Even though the (pC);c normalized sensitivity
coefficient at (y = 0) and L), are both negative and decrease with time, the k; CC normal-
ized sensitivity coefficients are different shapes with one decreasing (y = 0) and the
(y = Ly) values being positive. The difference in the k; cc and (pC):c normalized sensi-
tivity coefficients at (y = 0) and Ly is more pronounced as time increases and also after
22 seconds when the power to the heater is turned off. These sensitivity coefficients for
k; cc and (pC):c show that the experiment is well-designed because 1) the sensitivity

. e . . e .
coefficrents for ky, C are qurte different from those for (pC)cc and 2) the normalized

C

magnitudes are large (relative to the temperature rise of the experiment).

 

  
    

 

 

 

I I I l J _ l 1
A 2" / I \ (y = Ly) _
e / ”
a " I’ d
c: . / .
.2
U
E " .
—2~ -
g» - (y = 0) -
IE . .
...6— .-
. e e (y = 0)
. _ _ ky’ccaT/aky,“ .
. —— (0C)ECBT/B(0C)§c (y = Ly) .
_1G I I I I I

1 1
0 5 1 0 1 5 20 25 30 35 40

time (seconds)

Figure 2-7 Sensitivity coefficients for one-dimensional case 1023#508.2

36
2-5.0 Experimental Uncertainty
The estimated thermal properties are presented with confidence intervals indicating
the error in the estimates due to errors in temperature measurements. Only errors in the
temperature measurements are considered in such an analysis. Errors in other experimen-
tal measurements are not included in the confidence interval. Several other experimental
aspects influence the estimated thermal properties. These other experimental issues are not

the same nature as the errors in temperature. They are systematic errors and not random in

nature. The uncertainty in the estimate parameters (5), based on the uncertainty in the

experimental parameters, is calculated as

.. . A 1/2
~ ab 2 ab 2 3b 2
db = —d + —d + . . . —d 322
{(5.21) 0.. .2) +0., ..1 < >
where z = [z], z], . . . tsz are experimentally measured parameters. In the experiment to

estimate the thermal properties of the carbon-carbon, uncertainties in the measured heat
flux, location (depth) of the thermocouple, thickness of the carbon-carbon specimen,
effective thermal properties of the mica heater, and effective thermal properties of the
insulating material are considered. An experimental uncertainty in the thicknesses of the
mica heat and insulation are not included. This is because effective thermal properties
were experimentally estimated for these materials, implicit in the effective thermal proper-
ties is the prescribed thickness of the material. Consequently, uncertainty in the material
properties includes an uncertainty in the thickness. Experimental uncertainty for each
experimental parameter and its contribution to the total uncertainty in the estimated

parameters are given in Table 2-5.Uncertainty in the measured heat flux was computed as

37

Table 2-5 Experimental uncertainty for the two-dimensional experiments

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

P3328“ Uncertainty Contribution (g—Edz)
(pC):Cx10"° 14,“.
z dz
J/(m3°C) WI<m°C>
q 125 W/m2 0.033 0.049
L), 0.05 mm 0.016 0.027
y, 0.025 mm 0.014 0.015
kg,“ 20% 0.04 0.0
(pC)f,,,.ca 20% 0.06 0.06
kfm 20% 0.0 0.0
(pC)fm 20% 0.0 0.0
TOTAL 0.082 0.083
1
dq = {(gg-d/ijz + (3%;de +(%1d1j}2 (3-23)

where uncertainties in measurements of area (A), voltage (V), and current (1p) are
included.

There is not a dominant term in the uncertainty analysis. All experimental conditions
considered are of the same order. Overall, the uncertainty due to other experimental mea-
surements are excellent. The uncertainty represents a maximum of 2.4% in k; CC and
5.8% in (pC):C for the experiment at the lowest temperature and 1.6% and 2.8%, respec-
tively at the highest temperature. These values assume the uncertainties do not vary with

temperature, which is reasonable for these experiments.

Chapter 3

CARBON-CARBON THERMAL PROPERTIES: TWO-
DIMENSIONAL EXPERIMENTS

3-1.0 Introduction

In this chapter two—dimensional experiments are presented to estimate two compo-
nents of the thermal conductivity and the volumetric heat capacity for the carbon-carbon
composite. As in the previous chapter, effective thermal properties are estimated that rep-
resent the carbon matrix carbon fiber composite with a silicon-carbide protective coating.
The properties, k; CC, (pC):C, determined from the one-dimensional experiments are not
imposed on the two-dimensional solution. All the thermal properties are estimated simul-
taneously for the two-dimensional case. This permits a comparison between the results
from one-dimensional and two-dimensional experiments to demonstrate the consistency

of the methods.

3-l.1 Motivation

Interest in the solution of multi-dimensional inverse problems has gained momentum,
particularly in recent years. As the importance and wide-spread application of inverse
methods are realized, so too have the demands on the complexity of the problems that can
be solved. Such is the case for the application of inverse methods to the field of heat

transfer. Two examples are the estimation of the thermal properties of a material and the

38

39

determination of the heat flux at a boundary, both from experimental measurements. The
latter problem is the inverse heat conduction problem (IHCP), which has been the main
focus of research on multi-dimensional inverse problems in heat transfer. The application
of inverse methods to evaluate the IHCP or estimate thermal properties are closely related,

however.

A variety of methods which are used to solve the one-dimensional IHCP, have been
extended to the multi-dimensional case. Osman and Beck (1990), Hsu et al. (1992) and
Bass (1980) use methods based on the function specification method (Beck, et al., 1985).
Murio (1993b) has presented a mollified space-marching algorithm. The adjoint method is
employed by Jamy et al. (1991) and Truffart et al. (1993). Alifanov and Egorov (1985),
Alifanov and Kerov (1981), and Alifanov ( 1994) have presented formulations for iterative
regularization methods to solve the two-dimensional problem. The Monte-Carlo method
was investigated by Haji-Sheikh and Buckingham (1993). Solving the multi-dimensional

IHCP is further discussed in Chapter 4.

Although much energy has been focussed on the multi-dimensional IHCP, less work
has been afforded to the estimation of thermal properties using inverse methods for the
multi—dimensional case. Loh (1991) performed an experimental investigation for the esti-
mation of thermal properties, orthotropic thermal conductivity and isotropic volumetric
heat capacity, in a carbon-carbon composite. J arny et al. (1991) formulated the analysis to
estimate the thermal conductivity using a conjugate gradient method with an adjoint equa-

tion.

The lack of research on the multi-dimensional estimation of thermal properties may

be due to the small number of materials that display an appreciable anisotropy. Due to the

40

construction and advancement of composite materials, however, anisotropic thermal
properties are inherent in the composite; the magnitude of the anisotropy depends on the
type of materials. For the carbon-carbon material investigated by Loh the thermal
conductivity varied nearly an order of magnitude for directions normal and parallel to the
fiber direction. This anisotropic nature of the composite material requires multi-

dimensional inverse solution methods to accurately determine the thermal properties.

Although this chapter focuses on a laboratory method, the extension of a method to
the field, i.e. while the aircraft is on the runway or in the hanger, is of particular interest for
the carbon-carbon material because of the variability that the thermal properties demon-
strate. An in situ method also allows changes in the material properties to be tracked dur-
ing development and operation of the vehicle. The methods presented herein are not easily
extended experimentally to a field application, due to the practicality of instrumenting the
material. What is demonstrated, however, is the applicability of the analysis and algo-
rithms to determine the thermal properties given experimental data. More work on the

design and optimization of the experiments is required to move the methods to the field.

3-2.0 Experimental Aspects

The same experimental apparatus, shown in Figure 2-1, is used, except that only one
of three heaters is energized. The two-dimensional thermal model is shown in Figure 3-1.
The locations for the thermocouples are given in Table 2-1. The thermal model is mathe-
matically represented as three coupled problems. Although the techniques used herein can
be extended to temperature variable properties in a given experiment, the models given

below are for thermal properties that are temperature independent.

41

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11—
o I o o O L
n ceramic insulation " ms
W 6"
1 1
1 1,13 . . 2.14 : ll
/ composrte specrmen ( L
T/Ciifi 3,8 4,9 5,10 6,11 7,12 t y
’ 131 1:1 0 o o N A
i I" / , 77/ 1" <<<<<<<<<<(((((K/<<<<<<<<((((. '
11111100011011“ T x
q = 0
(10) L .
1/2 mica ””0"
heater
L assembly

 

Figure 3-1 Heat transfer model for two-dimensional experiments

The first problem models the mica heater/contact resistance

 

 

2 2
k8 a Tmica a Tmica _ C e Tmica _Lmica <)’< O 0
mica —2 +—2 — (P )mica a ,t>
33’ 33‘ 0<X<Lx

_ .- 91."... _ q(t)OSxSL1
Micaay )'=0 0 LI SX<LX
aTmica
E x=0L =0' ‘LmicaS-vso

TmicaU’ytO) = To’ -L aSySO .OstLx

mic

(3-1)

(3-28)

(3-2a)

(3-2b)

42

The second problem is for the carbon-carbon composite

2
ke 8 Ice
y,cc :;;§'

2
e ach_ Ce ach O<y<Ly O
”0517 —(p )“97 ’0<x<L ’t>
x
8T
—“ =0 OSySL
8x x=0L y

 

TQ(L%0)=7},O<y<LW,0<x<Lx

The third problem represents the insulation material

 

 

2 2
‘f Einuizim = (13C)? .33“ Li<Y<Ly+Li~s
"’3 Byz 3x2 ms at , O<X<Lx
aTins
0; _OL = 0,L),<y<Ly+L‘-ns
7L
g—ms = 0 OSxSLX
y y=L,-,,,+Ly

Tins(x,y,0) = To. Ly<y<Ly+LmS .0<x<Lx

The problems are coupled through interface conditions that assume perfect contact

k8

aTmrch, 0, t) _ ararx, 0:)

e

 

mica 8y - y, cc 0V

, armoc, L), z) _

ins

Tmica(x, 0, t) = T(C(x, 0, t)

_ e 0T(.C(x,L)_,t)
ay av

 

 

Tin5(x, Ly, t) = TCC(x, L), t)

,t>0

(3-3)

(3-4)

(3-5)

(3-6)

(3-7a)

(3-8b)

(3-8c)

(3-9)

(3-10)

(3-11)

(3-12)

43

Taktak (1992) investigated optimal experimental conditions to estimate the orthotro-
pic thermal conductivity for this geometry. However, his model specified a temperature
boundary condition at the backside of the composite, instead of (an approximate) insula-
tion condition which is used here. Taktak showed that for two sensors, the extremes
(x = 0, 7.62 cm) on the heated surface (y = 0) , was optimal for the condition that heat-
ing occurred over one-half of the surface. For this case, which is approximately insulated
at the backside of the specimen, includes the thermal effects of the heater assembly, and
heats over one-third of the surface, a similar outcome would be expected. Measurement of
the temperature on the surface where the heater is active and locations where it is not
active provides contrasting effects (sensitivity coefficients), which improves the accuracy

of the estimates. See the optimality criteria in Section 2-3.0.

3-3.0 Analysis Procedures

The same analysis procedures presented in Section 2-3.0 are used for the two-dimen-

sional case. However, now a two-dimensional heat conduction thermal model is solved.

Computer program PROP2D implements this inverse method to determine two com-
ponents of thermal conductivity and the volumetric heat capacity. The program was devel-
oped at Michigan State University by taking the finite element code TOPAZZD (Shapiro,
1986) and combining this direct problem solver with these parameter estimation methods,
to create a powerful algorithm. PROP2D allows for the estimation of the thermal proper-
ties for multiple materials, with possibly irregular geometries, from transient measure-
ments. The thermal conductivity can be orthotropic and temperature dependent thermal

properties are allowed.

3-4.0 Results and Discussion

A separate, independent set of one-dimensional experiments (Chapter 2) were per-
formed to determine the thermal properties of the mica heater and ceramic insulation in
the model (Figure 3—1). Effective properties were determined to account for imperfect con-
tact between the layers. Therefore, only the thermal properties of the carbon-carbon com-

posite are unknown in the model (Figure 3-1). Furthermore, one-dimensional experiments

e

were performed to determined the thermal conductivity normal to the fibers (ky, cc

) and
the volumetric heat capacity (pCzc) in the previous chapter. The one-dimensional results
provide initial estimates for the two-dimensional analysis and permit for a comparison to

demonstrate the accuracy and consistency of the methods; one-dimensional results are not

used as prior information in the analysis.

For the two-dimensional experiments, the analysis is more sensitive to the experimen-
tal conditions. For example, the magnitude and duration of the heat flux must produce ade-
quate response (sensitivity coefficients) for the sensors nearer the active heater, as well as
for the sensors farther from the active heater. This requires a longer heating duration than
that used for the one-dimensional experiments. The locations of the thermocouples must
also be known accurately, especially on the active heater where temperature gradients are
large. Although it is not too difficult to locate the position of the thermocouples in the car-
bon-carbon specimen, it is difficult to align the mica heater assembly relative to the speci-

mens since the heating elements are not visible in the heater assembly.

The two-dimensional thermal properties of the carbon-carbon composite determined
for temperatures up to 400°C are given next. Experiments were conducted at regular

intervals over this temperature range and analyzed assuming the thermal properties were

45

constant for the duration of an experiment, but varied between experiments. The measured
experimental data and details of the parameter estimation are presented and discussed for

a typical experiment.

3-4.1 Experimental Data

Experimental data for a test started at a temperature of 297°C (experiment case
1022$297) are shown in Figure 3-2. A sampling interval of 0.64 seconds is used to
acquired data for this experiment. The heating begins at approximately 16 seconds and
ends at 56 seconds. For experiments at lower initial temperatures the heating duration is
80 seconds. However, based on observation of the sensitivity coefficients and the criteria

for optimal experiments, “D-optimality”, (Beck and Arnold, Chapter 8, 1977) a shorter

 

 

 

 

 

 

 

 

325°G""'I‘1v1u..
sensors 3,8
320.04
4,9
(5 315.0-
E?
§ 510.0-
8.
“E’ 3050
[— o 1
300.0- 6?
:20000.0 §
295.0- Heat Flux ----- y = 0.9lcm _ ENE
\ — y = 0.0cm FIOOOOD é a
290'0""'I'1rj-..'o,o§
020"106080100120140 :15
time (see)

Figure 3-2 Experimental data for two-dimensional case 1022$297, see Table 2—1 for
sensor locations

46

duration, higher magnitude heat flux is selected to be closer to optimal experimental con-
ditions for determining the two components of the thermal conductivity and the volumet-
ric heat capacity. The complexity of this experiment is that two components of thermal
conductivity and volumetric heat capacity are simultaneously estimated. The experimental
conditions must be selected to provide information about all three effects. An alternative is
to conduct an series of experiments, each experiment providing information on one (or
more) particular effect. Then analyze the different experiments in a sequential manner.
Osman and Beck (1991) used such a procedure to estimate temperature dependent thermal
properties. For this model, a single experiment provides adequate information on all ther-

mal properties (effects) and a sequential procedure is not required.

The effect of the orthotropic thermal conductivity is seen by comparing in Figure 3-2
the temperature rise for the sensors at x = 3.81 cm (Figure 3-1) on the heated surface
(sensors 5,10) and x = 1.27 cm on the insulated surface (sensors 1,13). The larger thermal
conductivity in the x-direction results in a nearly instantaneous response at x = 3.81 cm
on the heated surface, while the sensor at x = 1.27 cm on the insulated surface has
approximately a four second time delay before responding. This delay exists even though
the sensors are approximately the same distance from the active heater (~10% difference,
0.835 and 0.914 cm from the sensor on the heated surface and insulated surface, respec-

tively).

Notice that temperature data are acquired after the heating ends. Continuing to
acquire data after stopping the heat flux results in better estimates. This is because it
causes the sensitivity coefficients to be of different shape from one another after heating.

These effects result in a more accurate estimation of multiple thermal properties based on

47
the criteria “D-optimality with constraints” (Beck and Arnold, pp. 459, 1977). Possible

heat losses in the experimental set-up can also be monitored with this data, although it
does not appear that there are significant losses in this experiment, since all temperature

SCIISOI'S Converge IO a constant.

3-4.2 Estimated Thermal Properties

Numerical issues are more important for the two-dimensional analysis than they are
for the one-dimensional. The mesh size and time step selected for the finite element
method can greatly influence the amount of time required to obtain a solution and the
accuracy of this solution. The mesh used for this analysis contains 525 (quadrilateral) ele-
ments. Twenty-five elements are along the 7.62 cm surface (x-direction) for all materials.
There is one element across the mica heater assembly and ten elements across each the
carbon-carbon specimen and ceramic insulation (y-direction). The computational time
step chosen was 0.64 seconds, the same as the experimental time step. A typical two-
dimensional analysis required 2-4 hours on a VAXstation 4000 Model 60 running VMS
V5.5—2, depending on the number of iterations and accuracy of the initial parameter esti-
mates; typically 4-5 iterations (16-20 direct solutions) were required. Although a detailed
investigation was not performed, the time step and mesh size were varied and shown to

have little influence on the resulting estimated thermal properties.

Two-dimensional thermal properties determined for the carbon—carbon composite are
summarized in Table 3-1. The experiment case number and initial temperature are given in
columns one and two. The next four columns present the estimated standard deviation (6)
and the two-dimensional thermal properties determined from the analysis with confidence

intervals. The final two columns give the duration and magnitude of the applied heat flux.

48

Table 3-1 Thermal properties estimated for the carbon-carbon composite from two-
dimensional experiments and analysis

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

EXP T1812: A k; cc k; cc (P C):cx 10-6
Case p 6 Heat Flux
(°C) (°C) W/mC W/mC J/m3c sec W/m2
1020&65 65 0.168 58.4 +/- 0.34 3.89 +/. 0.04 1.52 +7- 0.003 80 5276
10209111 111 0.183 61.0 +7- 0.39 4.17 +/. 0.04 1.74 +/- 0.004 80 6095
102119.157 158 0.163 60.7 +/- 0.34 4.24 +/. 0.04 1.93 +/. 0.003 80 6910
1021&205 205 0.210 61.8 +/. 0.33 4.55 +7- 0.04 2.13 +/. 0.004 80 8905
1021&256 256 0.202 61.6 +/. 0.30 4.73 +/. 0.05 2.34 +7. 0.004 80 8792
10228297 297 0.295 58.8 +/. 0.38 4.97 +/. 0.06 2.36 +/. 0.007 40 17304
10223352 353 0.277 58.7 +/. 0.37 5.07 +/. 0.06 2.56 +/— 0.007 40 17259
1022$403 403 0.256 57.4 +/- 0.34 5.09 +/- 0.05 2.66 +/- 0.007 40 17096

 

 

A plot of the thermal properties (ke kg and ( p C )2.) with estimated confidence

x, cc ’ y, cc ’
interval as a function of the initial temperature is shown in Figure 3-3. A F—test (Beck and
Arnold, Chapter 7, 1977) indicated a second order model (in temperature) for these prop-
erties. Because limited results are available, the relationship for k: cc remains linear. The

relationships determined with a least squares fit, are shown in Figure 3-3 and given by the

following relationships:

 

 

 

 

 

kg 387 193 LT‘ 0656 T_T' 2

= . . —. -13

e T— l
kt,“ = 60.8—2.07(T T J (314)

2_ 1
Ce 152 168 LT‘ 0543 FT‘ 2 10° 3 5
9-1-1 1—1 11x

C TZ—Tl TZ-‘Tl

where T1 = 65°C and T2 = 403°C. The thermal conductivity parallel to the fibers

49

 

 

 

 

 

 

8.0 . . . . - . - . . 8.0E+06
7.0 - e - 7.0E+06 _
. ky, cc .. 6
66.0- 2T1: =2: _._ ~6.0E+06§
é - = . E
E 5.0 — - 5.0E+06 .§‘
0
3? “ e ‘ 8.
E 4.0 _ kx, cc '1 4.0E+06 6
8 - « a
"O Q)
g 3.0 - - 3.0E+06 I
o . . .9.
"a 5

2.0 — - 2.0E+06
E . 9C; , g
..C: O
1‘ 1.0 - - 1.0E+06 >

000.0 5 100 ‘ 200 ‘ 300 ‘ 400 L 50(9-0E+O6

Temperature (C)

Figure 3-3 Temperature dependence of estimated thermal properties from two-
dimensional experiments

e
y, cc

(k.

x, cc) is 12 to 15 times as large as that normal to the fibers (k

). Testing was halted at
400°C due to the failure of the carbon-carbon specimen during subsequent one-dimen-

sional testing at higher temperatures.

In addition to estimating the thermal properties, PROP2D provides some means to
quantify the accuracy of the estimates. The estimated standard deviation, 6' computed in
equation (2-18), which is given in Table 3-1, provides an indication of how well the calcu-
lated temperatures match the experimentally measured temperatures. The magnitude of
the estimated standard deviation can be compared to the temperature rise of the experi-
ment, which is approximately 20 to 25°C. The magnitude of the estimated standard devi-

ation is within 1.2% of the maximum temperature rise for all the experiments and many

50

are less than 1%. The estimated standard deviation is largest for the last three experiments,

which applied a larger heat flux for a shortened duration.

There are other quantities that can be observed to demonstrate the accuracy of the

estimated properties. These quantities are the sequential estimates of the thermal proper-

ties, the temperature residuals, and sensitivity coefficients. These quantities are sensitive

indicators to provide insight to the estimation and experiment. Each is discussed below for

experiment 1022$297.

The sequential estimates demonstrate how the estimated properties vary as additional

experimental measurements are considered. Figure 3-4 shows the sequential estimates for

this experiment. A sequentially estimated property, at time ti, represents the outcome if

 

 

 

 

70.0 . l . . . . . . . . . 3.0E+06
A 50.0“ __ I‘M; ~2.5E+06
5; (90)..
\ 50.0-
E -2.0£+06
:5 40.0-
% -1.5E+06
g 30.0-
D
a -1.0E+06
g 20.0-
5

1 00- e '5.0E+05

—V¥ ky, CC
000 I T I I U 0.0E‘I'OO

 

 

 

time (sec)

I'jTl'l'l‘l
0 20 40 60 80100120140160

Figure 3-4 Sequential parameter estimates for two-dimensional case 1022$297

Volumetric Heat Capacity (J/m3C)

51

only data up to that time is used in the analysis and the results are linearized about the con-
verged parameter values. In other words, if the data is analyzed by adding one data set at
each time, it shows how the estimated properties change as one more data set is added to
the analysis. Initially the sequential estimates vary because there is not enough informa-
tion (data) to accurately determine the properties. However, as more data is considered,
the property estimates approach constants. If the experiment (or analysis) is ended at 80
seconds, the estimated properties would not differ significantly from the properties at 100
seconds. In general, for a good estimation the sequential estimates converge to a constant
and are steady with time. For this experiment the sequential estimates for k; cc, k; cc , and
(pC):C vary only 7.7, 3.0, and 1.6% over the final one-half of the experiment. The values
are quite large compared to the confidence intervals, which were 0.6%, 1.2% and 0.3%,

respectively. The confidence interval models error in temperature only. Uncertainty in

other experimental measurements account for the discrepancy.

The temperature residuals are related to 6 and are calculated as shown in equation (2-
18). They represent the difference between the measured and calculated temperature for a
particular time (ti) and sensor location (x j, yj) . The estimated standard deviation gives
an indication of the magnitude of the residuals; the sign and magnitude of the individual
residual can provide considerable insight. The residuals for this experiment (1022$297)
are shown in Figure 3-5a, b and c. Figure 3—5a presents the residuals for the sensors on the
active heater, Figure 3-Sb the residuals for the sensors on the heated surface but not on the
active heater, and Figure 3-50 the residuals for the sensors at the insulated surface. The
residuals for the sensors on the heated surface, but not on the active heater (Figure 3-5b),

are slightly correlated. The other residuals, Figure 3-5a and 3-50, are highly correlated and

52

 

Residual (C)
o
$

—0.504

 

 

----- Sensors 4,9 _

—1.00-

Sensors3,8
-1.50r1........ .

a 1 '
0 20 40 60 80 1 00 1 20 1 40

time (sec)

 

 

 

Figure 3-5a Temperature residuals for two-dimensional case 1022$297 on surface
below active heater

53

 

 

 

 

1.50 I | I I ' I ' l 1 T ' I ' I
1.00- _
,.. 0.50- -
8 . ' . .
-— ' 1 " I1 . , i ' \ . \ [Is I“: n d .’
g 0.00.. I 0* ‘ W '1 1 l r b 4
«1‘3 . 1 .
m I
.0501 _ Sensors 7,12 _
" _____ Sensors 6,11 ‘
_1‘OO- -— Sensors 5,10 3
. 4
“'1 .50 I I I I I I l I r I I I 1

O
N
O
‘5
O)
O
a:
O
—I
O
O
a...
N
O
E
.4.

time (sec)

Figure 3-5b Temperature residuals for two-dimensional case 1022$297 on the heated
surface, but not on the are of the active heater

54

 

    
 
 

 

 

 

 

1.00-
0.50-
g ., 'r
T; 0.00-
."2 ~ . ,
E I “‘0' Io‘,"1't.’.".‘t"0'-“o‘ ‘1‘ 0" 51v- ‘3’!
-0.50J
. ----- Sensors 2,14
-1.005
'4 Sensors 1,13
'1°50""'Iru-i...,
0 20 40 60 80 100 120 140

time (sec)

Figure 3-50 Temperature residuals for two-dimensional case 1022$297 at the
specimen/insulation interface

55

larger, 4—5% of the temperature rise on the heated surface and 8-10% of the temperature
rise at the insulated surface. It is not clear exactly why the residuals are so highly corre-
lated. Two possible reasons are the uncertainty in the location of the sensors with respect
to the heater and a non-uniform heat flux produced by the heater. First, there is some
uncertainty in the alignment of the heater assembly and the specimen. The heating ele-
ments are contained within an opaque fiberglass with a mica layer on the outside. Since
the fiberglass and mica extend beyond the heating elements, the alignment of the edge of
the heater (element) with the specimen is difficult. Second, the design of the heating ele-
ment has a gap between successive coils allowing for areas of localized heating. If the sen-
sors align with a gap, the actual heat flux will be larger than that calculated flux, which
assumes a uniform profile. Similarly, the sensors aligned with the heating element will
have a larger heat flux than the calculated uniform flux. The residuals for sensors away
from the active heater are less sensitive to the location and uniformity of the heat flux and

therefore are less affected.

To investigate the correlated residuals, two numerical experiments were conducted.
The first numerical experiment investigated an error in the location of the heater and the
second investigated a non-uniform heat flux; no other measurement errors were present.
These experiments were analyzed assuming no error in the location of the heater and a
uniform heat flux to investigate the effect on the residuals. Unfortunately, the results of
these numerical experiments were not conclusive. In particular, the two errors had
opposite effects on the residuals. An error in the location of the heater resulted in a
residual that had an opposite sign to the residuals for a non-uniform heat flux. The shapes

of the residual curves for the numerical experiments were quite different from the

56
2

experimental case (Figure 3-5), especially for sensors not at the location of the applied

heat flux.

As previously noted, observation of the sensitivity coefficients can provide insight to
the estimation problem. An advantage of parameter estimation, compared to other inverse
methods, such as gradient methods (J amy et al., 1991), is that the sensitivity coefficients
are computed in the analysis. Hence, no additional computation are required to observe
the coefficients. Observation of the sensitivity coefficients at this stage may be too late,
since the experiment is essentially designed and moving sensors or changing the heated
area is not easily done. However, some minor modifications may improve the accuracy,
such as changing the heating duration or magnitude. When possible, an analysis of the

sensitivity coefficients should be conducted prior to running the experiments.

Sensitivity coefficients provide insight to design experiments. In general, the normal-
ized sensitivity coefficients are desired to be large for each parameter and uncorrelated
(linearly independent) for different parameters. A sense of the magnitude of the sensitivity
coefficients is gained through normalizing the sensitivity coefficients. Normalization is
performed by multiplying by the parameters, resulting in unitsof temperature for the nor-

malized sensitivity coefficients. The normalized sensitivity coefficient for parameter 1] is

— 8T
X11 = n- (3-16)

A comparison is then permitted with the temperature rise of the experiment. For the cur-
rent experimental design, the normalized sensitivity coefficients are shown in Figure 3-6a-

c. Figure 3-6a shows the sensitivity to pCZC , ch, and Figure 3-6b and c shows the sensi-

tivity to the thermal conductivities k: c

.C and kin, X,“ and ka. Because the sensitivity

57

coefficients are normalized, a direct comparison of their magnitudes is possible. Some

observations are drawn from the sensitivity coefficient plots.

The most information is available on the active heater (sensor locations x = 0.89 cm
and 1.91 cm for y = 0). At these locations the sensitivity coefficients have the largest
magnitudes and therefore are most influential on the values of the estimated properties.

Notice that Xky undergoes a sign change across the specimen (in the y-direction) and ka

 

 

 

 

 

 

 

6 I r I I I I
4_ -
g 2- ,,,,,,, y=0.91cm ‘
B
C.
d)
'8 O -
35
O
0
U _2- ‘
3‘
IE
.2 _4_ -
d.)
m
-o
-6, 1
”<3
5 _3r _ t
Z -‘
‘10— sensor 3,8 I 1,13 2,14 d
—12 ' ' ‘ ‘ ‘ L
0 20 40 60 80 100 120 140

time (seconds)

Figure 3-6a Normalized sensitivity coefficient for volumetric heat capacity, X pC, for
two-dimensional case 1022$297

58

 

 

 

 

 

 

 

6 I . 1 I I I
- 1,13 _
4 . ......... y=0
A ....... =0.91cm
B 2- y _
a 712 2,14
0 0 ............... . ...... q
E
(3 6,11
3 '2" -
35 5,11
2 _4_ I
<1)
m
B
.5 '6“ _
E“
o -8— q
Z
sensor3,8
—10- -
-12 L l 1 1 1 1
0 20 40 60 80 100 120 140

time (seconds)

Figure 3-6b Normalized sensitivity coefficient for thermal conductivity in y-direction,
2,), for two-dimensional case 10225297

Normalized Sensitivity Coefficients (C)

-10

—12
0

Figure 3-6c Normalized sensitivity coefficient for thermal conductivity in x-direction

59

 

 

 

  

sensor 3,8

 
 

 

 

y 0
....... y - 0.9lcm .
20 40 60 80 100 120

time (seconds)

X k! for two-dimensional case 1022$297

140

60

undergoes a sign change along the specimen (in the x-direction). The implications of the
these sign changes are that there exist 1) a y-location within the body where the tempera-
ture is insensitivity to ky’ CC and more importantly 2) a x-location where the temperature is
insensitivity to k x, CC. The latter result is more important for this case where surfaces tem-
peratures are measured, because seemingly logical locations along the measurement sur-
faces y = 0, Ly may be insensitive to kx, CC. Although it is desirable to avoid locations
where the sensitivity coefficients changes sign because the sensitivity is quite small, sen-
sor locations that have sensitivity coefficients with opposite signs are a beneficial result.
This situation produces contrasting effects which improves the accuracy of the estimates.
Hence, the choice of the locations near the edges of the specimen (x = 0, 7.62 cm),
which have oppositely signed sensitivity coefficients, for the measurement locations. For—
tunately, the surfaces of the specimen (y = 0, Ly) are the most sensitive in that direction.
Finally, the experiment was halted after heating for approximately 40 seconds to make
comparable the magnitudes of the sensitivity coefficients for all thermal properties. In Fig-
ure 3-6b 70.) is approaching a steady-state value, while the sensitivity to X k! and XPC
have not. The sensitivity coefficient ch does not have a steady-state solution and will
continue to increase linearly with time. If the experiment heated for a longer duration Xk)‘
and ch would be much larger than X 1;,- Consequently, the estimates for the two former

properties would be considerably more accurate than the latter property.

3-4.3 Comparison of One- and Two-Dimensional Results
The two-dimensional results can be compared to the one-dimensional results in chap-
ter two for the thermal conductivity normal to the fiber and the volumetric heat capacity.

The results obtained for the one-dimensional analysis were used as initial estimates for the

61

two-dimensional analysis; however, the thermal properties were not constrained in the

two-dimensional analysis using the one-dimensional results.

The properties determined from the one- and two-dimensional experiments are com-
pared in Figure 3-7a and 3-7b. Overall, the results are quite close. The two-dimensional

and slightly lower for (pC):c. Recall that the pr0per-

C

results are slightly higher for k; C
ties were estimated assuming they were constant over the duration of an experiment and
therefore are applicable for a 20-25°C temperature range. Noting this, the one and two-
dimensional estimates for k; CC are within 6%, with the largest errors at the higher tem-
peratures where limited two-dimensional data is available. There is also a dip downward
for the estimates near the maximum temperature, an unexpected result that may indicate
the thermal conductivity is approaching a constant, but more testing (at higher tempera-
tures) is needed for verification. The results for (13C); have a maximum of 7% deviation

between the one and two-dimensional estimates.

The temperature dependence determined using both one- and two-dimensional

results, which were also indicated to be second order in temperature using a F-test, are

 

 

 

 

kg = 346+359 —208 3 ‘7
I I O -1
y'c (13 10) (7 3 10) ( )
Ce =ll4l+242 —08 '10 3 8
O O C -1
(P )CC (T3 7) 94( 3 ) x ( )

The relationship for k: CC has only two-dimensional results and is given by equation (3-
12). The estimated thermal properties, from both one- and two-dimensional experiments,
are within 4% and 6% of the relationships given in equation (3-17) and (3-18) for k: CC

and (pC):C.

62

 

 

2 O— . 2DExp
O lDExp

Thermal Conductivity (W/mC)

 

 

 

' ' 1 1 1 1 1 ' 1
O 1 00 200 300 400 500
Temperature (C)
e
11) (PC )Cc

 

3.0E+06—
2.5E+06-
2.0E+O6-—

1.5E+06~

 

1.0E+06— 0 21) EXP -

« O 11) Exp -
5.0E+05— -

Volumetric Heat Capacity (W/m3C)

 

 

0.0E+00 v 1 1 1
0 100 200

 

I I I I

1 1 1 1
300 400 500 600

Temperature (C)

Figure 3-7 Comparison of thermal properties estimated from one- and two-

dimensional experiments a) thermal conductivity k: CC and b) volumetric specific heat

(0C);

63

Table 3-2 Experimental uncertainty for the two-dimensional experiments

 

 

 

 

 

 

 

 

 

 

 

 

 

Paigeter Uncertainty Contribution (gig-dz)

(pC1:.x10“° k2... k5,...
z dz
J/(m3C) W/(mC) W/(mC)

q 3%*q W/m2 0.086 1.9 0.15
Ly 0.05 mm 0.014 0.43 0.066
y1 0.025 mm 5.5 13-04 0.14 0.016
x J. (i=1,2. J) 1.0 mm 8.0 13—04 0.72 0.20
kg,“ 20% 0.018 0.55 0.073
(003,“, 20% 0.020 1.2 0.14
kfm 20% 0.0 0.0 0.0
(0C); 20% 0.0 0.0 0.0
TOTAL 0.091 2.5 0.30

 

 

 

 

 

 

3-5.0 Experimental Uncertainty

Using the techniques described in Section 4-5.0 for the one-dimensional properties,
the experimental uncertainty for the two-dimensional thermal properties is assessed. For
the two-dimensional experiment uncertainties in the measured surface heat flux, thickness
of the carbon-carbon specimen, y-location and x-locations of sensors on the surface near-
est the heater (experiment is insensitive to locations specified for sensors on back surface
of specimen), thermal properties of mica heater, and thermal properties of insulation are
considered. As in the one-dimensional analysis, the uncertainty in the thicknesses of mica
heater and ceramic insulation are not considered because they are implicit in the effective

properties estimated for these materials. Table 3-2 gives the uncertainties in the estimated

64

thermal properties due to the uncertainties in the measured quantities. For the two-dimen-
sional case there are three important uncertainties. The uncertainty in the measured heat
flux, x—location of the sensors, and the thermal properties of the mica heater are dominant
terms. All terms contribute significantly to the overall uncertainty, but vary depending on
the thermal property. The heat flux is important because there is a large uncertainty in the
area associated with locating the heater relative to the sensors. The x-location of the sen-
sors are important due to alignment issues. The two-dimensional uncertainties are larger
than one-dimensional results. However, uncertainties are still a maximum of 6.0, 4.2, and

7.7% of(pC)" k" and kf’

CC , x, CC, v, CC, respectively.

Chapter 4

SOLUTION OF IHCP USING A GRADIENT METHOD WITH
ADJ OIN T EQUATION APPROACH

4-1.0 Introduction

Estimating the conditions at the surface of a conducting body from internal measure-
ments is typically called the inverse heat conduction problem (IHCP). Inverse describes
this type of conduction problem because conditions at the boundary or surface of a body
are estimated using internal measurements. Whereas a direct conduction problem, some-
times called just the direct problem, uses conditions specified on the boundary to compute
the internal temperature. While the direct problem is generally a well-posed problem, the
inverse problem tends to be ill—posed and very sensitive to measurement errors. A sequen-
tial method to solve the IHCP is discussed in this chapter.

Applications of the IHCP are abundant. A classic problem is determining the heat flux
at the surface of space vehicles during re-entry. A related problem is the estimation of the
heat flux at the leading edges of hypersonic vehicles. Other examples are estimating the
heat flux, or convection coefficients during quenching, which enables control in the result-
ing material properties. During casting processes, estimates of the heat flux can be used to
the control the cooling rate and movement of the cooling front to minimize residual

stresses and void formation. Measurements on the outside of the cylinder wall of a internal

65

66

combustion engine can estimate the heat flux in the combustion chamber. Additional
examples are cited in Kurpisz and Nowak (1995).

Several methods are applied to solve the IHCP. Function specification (Beck et al.,
1985), Tikhonov regularization (Tikhonov and Arsenin, 1977), gradient methods (Ali-
fanov, 1994 and Ozisik, 1993), and mollification (Murio, 1993) are more frequently cited
methods. However, other approaches also applied are, dynamic programming (Busby and
Truijillo, 1985), Kalman filter (Tuan et al., 1996), Monte Carlo method (Haji-Sheikh and
Buckingham, 1993). In addition, combining methods is useful. Beck and Murio (1986)
formulated a combined function specification and Tikhonov regularization method, simi-
lar concepts were used in Osman et al. (1997). Jamy et al. (1991) combined an iterative
search method with Tikhonov regularization. Since this dissertation focuses on multi-
dimensional problems, the literature review will narrow its attention to the two-dimen-
sional problem. Refer to books by Alifanov (1994), Beck et a1. (1985), Hensel (1991), and
Kurpisz and Nowak (1995), for a comprehensive survey of the literature.

Many methods applied to the one-dimensional problem have been extended for the
two-dimensional problem. Function specification and gradient methods have received the
most attention. Function specification specifies a functional form for the heat flux over a
future interval (Beck et al., 1985). Specifying a functional form over the future interval
provides regularization to stabilize the ill-conditioned problem (Lamm, 1995). In conjunc-
tion with specifying a function form, function specification solves the problem in a
sequential manner.

Several researchers have investigated the two-dimensional application of the function

specification method. It was applied to estimate spatially and time varying convective heat

67

transfer coefficients, Osman and Beck (1989, 1990), and surface heat flux, Osman et a1.
(1997). A boundary element method was coupled with function specification to investigate
multi-dimensional problems by Zabaras and Liu (1988). Hsu et al. (1992) applied a finite
element method to solve the general two-dimensional problem with inverse methods simi-
lar to function specification.

Gradient methods, which typically apply a conjugate gradient iterative scheme, cou-
ple iterative or Tikhonov regularization to stabilize the solution and solve the multi-dimen-
sional problem. Iterative regularization depends on the slowness or “viscosity” of the
solution and uses the iteration index as the regularization parameter. Several papers use
iterative regularization; see Alifanov and Kerov (1981), Kerov (1983), and Alifanov and
Egorov (1985). Additional investigations using gradient methods, but not iterative regular-
ization, are given in Zabaras and Yang (1996), Reinhardt (1996a.b) and Jamy et al. (1991).

Others have studied the multi-dimensional inverse problem with a variety of
approaches. Tuan et al. (1996) uses Kalman filter to develop an on-line algorithm. A new
method, called “direct sensitivity coefficient,” is claimed by Tseng and Zhao (1996) and
Tseng et a1. (1996). Pasquetti and Le Niliot (1991) employ Tikhonov regularization with
the boundary element method. An adjoint equation approach to compute sensitivities and
relate measured temperature to unknown surface conditions is used by Hensel and Hills
(1989). A Monte-Carlo method is given by Haji-Sheikh and Buckingham (1993). Mollifi-
cation with a space marching technique is used by Murio (1993) and Guo and Murio
(1991). Busby and Trujillo (1985) use dynamic programming. An analytical solution is

developed by Mosaad ( 1995) and transform methods are used by Imber (1974, 1975).

68

Function specification and gradient methods are popular approaches. An advantage of
the function specification method is that the problem retains the causal nature, represented
with a Volterra operator (Lamm, 1995), and can be solved sequentially with possible sav-
ings in computational time and memory. The method, however, requires assuming a priori
information about the (unknown) function. Gradient methods require no a priori informa-
tion but typically solve for a function on the whole time domain, not taking advantage of
the causal nature of the problem.

The demonstrated success of gradient methods and the efficiency of a sequential
implementation suggest a sequential implementation of a gradient method would be a
powerful combination. A method that sequentially implements a gradient scheme, using
an adjoint equation approach, is proposed and to be developed; additional stability is intro-
duced by including Tikhonov regularization. The method is anticipated to benefit from not
requiring a prescribed functional form, which is particularly useful on space, and the com-
putational aspects of a sequential implementation. Other researchers who have proposed a
sequential implementation are Reinhardt and Hao (1996a,b) and Artyukhin and
Gedzhadze (1994); however, past implementations have addressed the one-dimensional
problem with a very limited investigation of the method. This dissertation addresses the
one- and two-dimensional problems for the sequential implementation. Although applying
the sequential gradient method to linear one-dimensional problems is probably not needed
or recommended, addressing the one-dimensional problem provides valuable insight to
the method. This is the first known two-dimensional implementation of a sequential gradi-
ent method and the most comprehensive, to the best of the author’s knowledge, study of

the gradient methods.

69

Derivation of the equations to solve the IHCP using a conjugate gradient search
method with an adjoint equation approach are presented in this chapter. Gradient methods
that apply an adjoint equation approach have been widely discussed in the literature. In
particular, the Russian community has championed these methods (Alifanov, 1994). How-
ever, the mathematical basis of the method seems to have limited application of the meth-
ods (at least for engineers), more so in the USA than elsewhere. Reasons for the limited
use of the method in the USA can be partially attributed to the strong advocacy of the FSM
there (Beck et al., 1985). However, more recently the power of the methods, particularly
for multi-dimensional problems and problems with many estimated components in the
author’s opinion, has been realized. In particular, Ozisik (1993), is an advocate of gradient
methods. In addition, the methods have been presented in a less mathematically rigorous
form (Jamy et al., 1991, Lamm, 1990, and Jamy and Beck, 1995).

Section 4—2.0 to 4-6.0 develop the describing equations to solve the general multi-
dimensional IHCP. Equations are developed in differential form and presented in the form
of an algorithm in Section 4-7.1 for a whole domain solution. A sequential implementa-
tion is discussed and an algorithm is presented in Section 4-7.2. Derivations up to this

point in the chapter assume that the unknown function resides in a prescribed function

space E, which is an infinite-dimensional space. Restricting the solution to a finite-dimen-

sional space, Rn , is addressed in Section 4-8.0. The chapter is concluded with a section to

summarize the chapter and provide some insight to the method.

70
4-2.0 Problem Statement of the Multi-Dimensional IHCP

A schematic of the general multi-dimensional IHCP is shown in Figure 4-1. In the
analysis that follows, the conduction problem is assumed linear, i.e., thermal properties do
not depend on temperature. Extending the gradient methods for non-linear problems is
discussed in Artyukhin (1996) and Loulou et al. (1996). For the sequential implementation
to be developed, it is possible to temporarily linearize the problem and not consider it non-
linear. Assuming the thermal properties are independent of temperature, the problem is

mathematically formulated as follows:

V. (k(r)VT(r, t)) = pC(r)[—aa—IT(r, t) + v - VT(r, 0], at)": :nSgt:) (4-2.1)
_k, 97170", t) + h, T(r, t) = fi(r, t), (r) on 1“,, (i =1’ 2’ 3) (42.23)

(t0<tStf)

(T=f](r7 1))

     
  

Y2(t) Yj(t)

   

(qu2(rrt))

(610,1): (0

Figure 4-1 Schematic of multi-dimensional general IHCP

71

_k(r)-—87T(r, t) = q(r,t), (r) on F4 (4-2.2b)
3" (tO < t s tf)

T(r, t0) = Tom (4-2.20)

where k(r) and pC(r) are the thermal conductivity and volumetric heat capacity. The
spatial domain 9 is moving at a constant velocity of v and the symbols

1“,- (i = l, 2, 3, 4) representing the domain boundaries. The outward pointing normal
vector is denoted fr. Functions f1(r, t), f2(r, t), f3(r,t), and T0(r) are assumed
known. Boundary coefficients kl. and hi are specified to form different boundary condi-
tions, i.e. k1 = 0 and h1 = l specifies a temperature boundary condition (first kind),
k2 = k(r) and h2 = 0 a flux condition (second) kind, and k3 = k(r) and h3 = h(r)
a convective condition (third kind). Surface F4 has the unknown heat flux, q(r, t) , to esti-

mate. Excess information of transient temperature measurements exists within the body

(or at the surface of the body) at locations (r =dj) j = 1, ....... J . These measured tem-

peratures are denoted
Y(dj,tk) = YJ-(tk) (42.3)

and are available at discrete times tk, 0 < tk S 1f, k = l, 2 N. The analysis is first devel-

oped for measurements being at discrete locations but continuous on time.
The objective is to estimate q( r, I) using the measured temperatures. Hence, the heat
flux q(r, t) is estimated such that the calculated temperatures match the measured tem-

peratures, T(dj, tk;q) = Yj(tk) , where T(dj, tk;q) is the solution of equation (4-2.l) and

72

(4-2.2). Due to the ill-conditioned nature of this problem, the matching is accomplished in

a least squares sense by minimizing the function

J
%E £’1T(d,-,t.;q)— 13-10124: + gaTfi’jnrqu, t) _qpn.(,, 012.1, d,

J(q) = j=1

 

f

P J

v

15 J R
(4—2.4)

The last term in equation (4-2.4), J R, serves as regularization to stabilize the ill-condi-

tioned problem; the function q [NI-(s, t) is known and could be zero. Regularization of this
form is similar to zeroth order Tikhonov regularization (Tikhonov and Arsenin, 1977).
Schemes to minimize J (q) in equation (4-2.4), which use iterative methods such as
steepest descent and conjugate gradient, require the gradient of J (q) . Methods to compute
this gradient depend on the function space where q(r, t) is assumed to reside. Two possi-

bilities are a finite dimensional space and an infinite dimensional space. For the infinite

dimensional problem a priori information is not required concerning the (unknown) func-
tion q(r, t) . However, computation of the gradient requires solving two additional prob-

lems, which are the adjoint and sensitivity problems. For the special case when a priori
information is available (or assumed) concerning the function, the problem is considered
finite dimensional and standard differential calculus is used to compute the gradient.

The approach considered here is the more general infinite dimensional problem. The

heat flux is assumed to be in the function space E on F4 x (to, if) with the scalar product

defined by

(Z,,2,)E = I? [2,0, I)Zz(r, t)dr d: (42.5)
”F,

73

for functions Z1(r, t) and Z2(r, t). The associated norm or function Z (r, t) is

HZ“: = JJIJIZU, r)|2dr dt (4-2.6)
to 1“,,
One possible choice for the function space is E = L2 , all square integrable functions on

F4 x (0, if) . The finite dimensional problem, i.e. E = R" , is discussed in section 4-8.0.

4-3.0 Gradient Calculation

The gradient of the function J (q) in equation (4-2.4) is required for gradient search
methods. In all derivations that follow it is assumed that the function J (q) is given by
equation (4-2.4). This is an important point because J (q) is quadratic, since temperature
is a linear function of the heat flux q. Application of the methods for non-linear problems
is discussed in Artyukhin (1996) and Loulou et al. (1996). Because a functional heat flux
is assumed, (q e E), the mathematics are more complex. The gradient of J (q) at q,

denoted VJ (q) , is related to the variation of J at this point by the general relationship

J(q + Aq) — J(q) = (VJ(q), Aq)E + (terms nonlinear in llAqIIE) (4-3.1)
Although more mathematics are required to correctly interpret this gradient, it is analo-

gous to the generalization of a Taylor series expansion in standard differential calculus.
For the function J(q), the Gateaux differential with increment Aq, henceforth
referred to as the directional derivative at q in the direction Aq, denoted D AqJ (q) , is

defined as (Luenberger, 1973)

DAqJ(q) = Ji-TOJ(q+uA:)—J(q) (4.32)

The directional derivative D A q] (q) is therefore related to the gradient

74
DAqJ(q) = (VJ(q), Aq)E (43.3)

The directional derivative is computed using equation (4-3.2); however, identifying the
gradient requires the derivation of two additional problems, typically called the sensitivity

and adjoint problems.

To compute the directional derivative of J (q) the limit in equation (4-3.2) is evalu-

ated
J(q + qu) -J(q) =

J
l I . 2 l f 2
5.21J:0 [T(dj, t,q + qu) - Yj(t)] dt + 2a7£0Jr.[q(r’ t) + 11440“, t) - qpn-(f, t)l dr 0’!
J:

J
1 2 1 2

-§ 2 flfldj, t, :61) - Yj(t)] dt— 2aT£fIr [q(r, t) —qpn.(r, t)] dr dt(4-3.4)
jzl 0 0 4

After expanding the terms, dividing by u, and some manipulation, equation (4-3.4) is

rearranged as

J(q+qu) -J(q) :
u

 

.1
$2 £:{1T<d,.t;q+qu)+T<d,.t;q>—2Y,-0>1 u

[T(dj, m + 114(1) — T(d,» (1)1}dt
j=l

 

2
2 9 — ' ’ A 1 A ’
“Lialeflr { [610 I) qp,,(r 0111 610 0+1” q(r 1)] }dr W435)

u

The resulting change in the temperature caused by a variation of magnitude ItAq in the
heat flux is related to the directional derivative of the temperature. The directional deriva-

tive of T at q in the direction Aq, evaluated at (r, t) is

 

DAan, t;q) = lim T(r’ “‘1 + “A” ‘ T(r, ”1) (43.6)

u-—)0 [1

Taking the limit of equation (4-3.5) as It —> 0 while using the definition of the directional

derivative of temperature in equation (4-3.6) gives

75

J
DAquq> = 2 J2me. t;q) - Y ,(r)1DA,T(d,-. t;q)dr
j: 1

+ aTJ’fI [q(r, t) — qpnrr, t)]Aq(r, t)dr (1:14-37)
to 1‘4
where the left side of equation follows from the definition in equation (4-3.2). This equa-

tion shows that the directional derivative of the J (q) is a function of the directional deriv-

ative of the temperature D A qT(d j, t;q) . The directional derivative of temperature, i.e., the

variation in temperature due to a variation in the heat flux, is typically referred to as the
sensitivity function and defined as
0(r, t;Aq) E DAqT(r, t;q) (4-3.8)

The sensitivity function 0(r, t;Aq) represents the variation in temperature T(r, t;q) due

to the presence of an input heat flux function Aq.

4-4.0 Sensitivity Equations

The describing problem for the sensitivity function is derived from the definition of

the directional derivative in equation (4-3.6). Evaluating the direct problem in equation (4-
2.1 and 4-2.2) at q + ItAq gives

V ' (k(r)VT(r, I:q+11461))

= pC(r)[§—T(r, t;q + ItAq) + v - VT(r, t;q + ItAq)], (r) in Q (4-3.9)
at (to < t S If)

(r) on 1",, (i =1, 2, 3)
(t0<tStf)

-k,- Egg-T(r, t;q + qu) + h,- T(r, t;q + qu) = fi(r, t), (4-3.10a)

(r) on F4

(tO<tStf)

—k(r)§fi-T(r, t;q + qu) = q(r, t) + qu(r, t), (43.1015)

76
T(r, 0;q + ItAq) = T0(r), (r) in Q (4-3.10c)
Defining Tp = T(r, t;q + qu) in equation (4—3.9) and (4-3.10), subtracting the direct

problem in equation (4-2.1) and (4-2.2) with T = T(r, t;q) , and dividing the result by It

 

 

 

 

 

 

 

 

gives
T —T T —T T —T '
V-(k(r)V( D = pC(r)[3( 1‘ )+y-V( 1‘ j] , (”mg (4311)
at 11 ll (t0<tStf)
T —T T —T , -_._.
—k,. a‘( “ )+hi( ’1 )z 0, (”on IT“ 1’2’3) (4-3.12a)
3" ll 11 (t0<tStf)
T —T
-k(r)-§—.-( u )= Aq(r,t). (”on F4 (4-3.12b)
3" 11 (IO<tStf)
T —T
( “ )l = 0,(r) inn (4-3.12c)
11 (:0

Taking the limit as It -—> 0 of equation (4-3.l 1) and (4-3.12) while using the definition in
equation (4-3.6) and the definition of the sensitivity function 0(r, t;Aq) E DAqT(r, t;q) ,

the describing problem for sensitivity function is

V-(k(r)V0(r,t)) = pC([)[—a—0(r,t)+v - V0(r,t)], (r) i“ 9 (43.13)
at (t0<tStf)
_kla__0(r, t)+h 9( )= 0’ (r) on 1“,,(i=1,2,3) (4-3.143)
(t0<tStf)
_k(r)— "0(r, I): Aq(r,t), (r) on F4 (43140)
(t0<tStf)

0(r, 0): 0, (r) in Q (4-3.l4c)

77

The describing equations for the sensitivity function have a similar form to the equations
derived for direct problem. The differential equation and boundary conditions have the
identical form. Differences between the direct and sensitivity problems are that the sensi-
tivity problem has homogeneous known boundary conditions, equation (4-3.14a), it has a

zero initial condition equation (4-3.14c), and its driving term is the variation in the heat

flux Aq(r, t) on surface F4, equation (4-3.14b). Similarities in the form of the two prob-

lems (direct and sensitivity) allow for savings in the calculations required for the numeri-
cal solution of these two problems. It is shown in the subsequent section that the

describing equations for the adjoint function have a similar form.

4-5.0 Adjoint Equations

To obtain the derivative, VJ (q), an adjoint problem is derived which will allow the

directional derivative to be written as the scalar product
DAqJ(q) = (VJ, Aq)E (44.1)
The current expression for D A q] (q) in equation (4-3.7) is not separable into a scalar prod-

uct that contains Aq, i.e., (function, Aq) E where terms that comprise “function” would

then be the gradient. Hence by developing an adjoint problem the directional derivative is
written as a scalar product and the gradient is identified.

To derive such an adjoint problem the method of Lagrange multipliers is used. As
posed, the problem is to determine a solution q*(r, t) that minimize J(q), subject to
q*(r, t) satisfying the direct problem, equation (4-2.1) and (4-2.2). Equations (4-2.l) and

(4-2-2) are rewritten such that A[T(r, t;q)] = 0, where A is an appropriate operator. A

78

necessary condition for q*(r, t) to satisfy these conditions is that there exist Lagrange

multipliers, Mr, I) , such that the Lagrangian function
110., T, q) = M. q) + 1140-, 041m. t;q)ldr dr (44.2)
i Q

is stationary at the solution, q*(r, t). This implies that the directional derivative of
A0», T, q) is zero at q*(r, t). (In equation (4-4.2), the integrations must be interpreted in
the correct sense for the boundary conditions associated with operator A .)

Using Lagrange multipliers permits the temperature T(r, t;q) to be considered inde-

pendent of the heat flux, q(r, t) . This is possible because the interaction between the tem-
perature and heat flux is forced in the constraint, eliminating the need to consider the
temperature as a function of q.

Define a Lagrangian function which includes the criterion J (q) and Lagrange multi-

pliers (w, 11,-, AU) associated with the imposed constraints that the heat flux must satisfy

the direct problem gives
f 0T

4011.11,. k... 7.4) = J(T.q>+ f jv(r.t){V - (k(r)VT)—pC(r)[3-t +v - VT]}dr 4:
[On

3
+ 2 Eryn”, ”(19. %%_hiT+fi(r, t)}lr dt
”r

I31 ,‘

+ I? In4(r, t)(k(r)g% + q(r, 1))dr dt + I710(r)[T(r, 0) _ Tom] d, (4.43)
01-44 Q

Applying Green’s theorem and integrating by parts on time, the first integral in equation

(4-4.3) can be rewritten as

79

JII‘IWO', t){V . (k(r)VT) - pC(r)[%lt~ + v . VT]}dr dt =
0::

£’I{V . (k(r)V\V) + [pC(r)g—}V+ v - V[pC(T)Wl]}T 61" d’
052

+ ffl‘l’kUg—g— “HT—313w dt-LIIPCUNITU . fildr dt
0F 0F
‘jPCWT )1?= tow-(4-4.4)
n

where F = F1 u F2 u F3 U F4. Substituting this integral into the Lagrangian in equation

(4-4.3) gives

4111.11.40.14) = J(T.q>+f’j{V~1k<r)W)+ 960) 31,-" +v - Vipcrr) w1}r dr 4:
”Q
I 3T 3 f .. I
+ “[1405;— 1.6113434, .11 _ 1‘0 ijcvwnv - n>dr dt— #ch >1;= ,Odr

+ 2]: 111,0: )k,[g% h T+f.,(r t)]dr dt+J:: [1140, t)[k(r)gI-fi +q(r 0].), 4,

i=1

+ I;o(r)[T(r, 0) - T0(r)] dr (4—45)
9

The Lagrangian has now been rewritten such that the differential operators in the first inte-
gral apply to the Lagrange multiplier 01, instead of temperature T. The directional deriva-

tive of the Lagrangian function A, for fixed multipliers, is

SAW, 11;" A.» T, q) = DATAW, 111, AG, T, q) + DAqu, ‘11, 10, T, q) (4-4.6)

where D ATA( w, 11,-, 10, T, q) is the directional derivative at T in the direction AT
evaluated at (11!, n i, 21.0, q) and D A qA(w, m, 3.0, T, q) is the directional derivative at q in
the direction Aq evaluated at (111, T1,, 71.0, T). An analogous expression, similar to equation

(4-4.6), can be written for the directional derivative of J(T,q). These directional

80

derivatives can be obtained by using a similar methodology to that given in equation (4-
3.2).

The directional derivative of the Lagrangian criterion for fixed multipliers is

SAW, 11,-, to, T. q) = 51(T. q) + Jj’j{V - (k(r)V\v) + pC(r)? + v - leC(r)w1}AT dr dt
on
B 8 . t
+ [Thane—fin — k(r)AT3—¥]dr dt — £’IpC(r)wAT(v - n)dr dt — ({pCWAT )1]= I0dr

+ 2f: [11,0 t)l:k-aan —AT+— —h.AT]dr 811+]: [1140, t)[k(r)aan —AT,+Aq(r t)]dr d,

i=1

+ jl,(r)1AT(r, 0)] dr (44.7)

9
The solution of interest is for the case that the temperature satisfies the direct problem, i.e.,
T = T(q) . Hence, for this case AT 2 0 in equation (4-4.7), therefore the boundary con-
ditions on sensitivity problem, equation (4-3.l 1) and (4-3.12), can be used to simplify this

equation. Introducing the boundary and initial conditions from the sensitivity problem, the

Lagrangian criterion in equation (4-4.7) simplifies to

6401!. T14). q) = 611T<q>.4>+f’1{V- (k(r)Vw)+ 96013;" + v - V10C<r>w1}e dr 4:
()Q

3

2130— —JJ I:I{1I/Aq(r, ’)+[k(5)a-7 +pC(r)(v n)\|l]0 }dr dt

i=1

-ij(r)(\1I8 )1, dr (4-4.8)
Q f

where,

 

(i = 1) 13C! = K’Hkm 32—pC(r)(v-n)0]wdr dt (449)
(Ir

1

81
(i¢ 1) IBC, = JTJI— k(r)%%’+[hi—pC(r)(v~fz)]w]6dr dt (4-4.10)
()ri

Notice that the criterion is no longer a function of the multipliers “i and lo since the

terms for these multipliers are zero from the boundary conditions on the sensitivity prob-

lem. It can be shown that

51(T(q), q) = DAqJ(q) (4-4.ll)

Defining the residual for the difference between the measured and calculated temperatures
(ej(t) = YjU) - T(d,» I, 4)) in the directional derivative DAqJ(q) , equation (4-3.7), the
directional derivative of the Lagrangian is

SAW, T(q), q) =

J
_ 2 £II[ej(t)5(r—dj)]6 dv dz + an’jr [q(r, z) —qp,,.(r, t)]Aq(r, z) dr dz
j=l 09 t0 4

+ij [{V - (k(r)Vw) + pC(r)-git" + v ~ V[pC(r)w]}9 dr dt
"Q
3

+ Z IBC, —J:f I{qu(r, t) + [k(r)g—:i: + pC(r)(v - ii)\|l]6} dr (1!
or4

i=1

-IPC(f)(W9 )|t=t dr(4-4.12)
Q I

Rearranging the terms in equation (4-4.12) to group similar integrals gives

8A(W’ T(q)? q) =
J

K‘HV ' (k(rw‘”) + 9C“)? + V ' VIPC(r)vl - Zlej(z)6(r_dj)}e dr dz
J:

+ GTE-(IL [q(r, t) — qpn-(r, t)]Aq(r, t) dr (1'!

3
+ 2 13C, ‘ij j{qu(r, z) + [k(r)%%’+ pC(r)(v - fi)w]9} dr dz
“F4

i=1

+ JPCUXWG )|t=t dr (44.13)
9 j

82

Fix the Lagrange multiplier, w(r, t), as follows:

J
V . [k(r)Vw]+pC(r)%V +v . VIPC(r)\Vl— z e,(t)8(r—d,-> = n.
j=1

(r) in 9
(44.14)
(r0 < r s tf)

-kz—a,-wtr,z>+Ihz—pC(r><v-fi>1ww) = 0. (')°"F"(’=1’2’3) (44.15»
3" (t0<tStf)

k(r) 5%“!0, z) +pC(r)(v - it) Mr, z) = 0, (r) on F4 (4-4.15b)
(to < t s tf)

w(r,t)|1_l = O, (r) in Q (4-4.15c)
" I
By defining the Lagrange multiplier as done in equation (4-4.l4) and (4-4.15), the

Lagrangian function in equation (4-4.l3) reduces to

SAW, T(q),q) = -£fJWAq(r, I) dr 61”“sz [q(r, t)-—q,m-(r, t)lAq(r, 0dr dt
0r4 0 4
(44.16)

Instead of fixing the Lagrange multiplier as prescribed in equation (4-4.14) and (44.15),
the negative of the Lagrange multiplier is calculated. Since the Lagrange multiplier is
homogeneous except for the residual term, changing the sign of the residual term only
changes the sign of the multiplier. Define the adjoint function as this Lagrange multiplier,
with describing equations

J
V . [k(r)V\p] +pC(r)%V+v - V[pC(r)w] + 2 ej(t)8(r—dj) = O,
j=l
(r) in Q
U0<t3rfl

(4-4.l7)

83

-kz 13W, t)+[hz-pC(r)(v-ii)] Mr. I) = 0. (r) on no =1’2’3) (4—4.18a)
3" (t0<tStf)
k(r) BAH", t) + pC(r)(v ~13) w(r, t) = O, (r) on F4 (4-4.18b)
3" (t0<tStf)
WU, t)|,=,f = 0. (r) in Q (4-4.18c)

By defining the Lagrange multiplier in this manner, the directional derivative of the

Lagrangian reduces to

MW, T(q), q) = 1:: r{1113610, t) dr d! + arfiflrfqm t) —quz(r, t)]Aq(r, t)dr dt

(4-4. 19)

To determine the relationship between the directional derivative of the Lagrangian
and the function J, consider that temperature T satisfies the equations given in equation

(4-2.1) and (4-2.2), denoted T(q) , and the Lagrangian criterion simplifies to

AW, T(q), q) = J(T(q), q) (4-4.20)

This relationship is seen from equation (4-4.5) where all integral terms are zero for

T = T(q) . Assuming fixed multipliers, the directional derivative of equation (4-4.20) is
5A(\V, T(q),q) = 51(T(q),q) = DAqJ(q) (4-421)

Since equation (4-4.19) is also satisfied for T(q) , it can be rewritten using the result of

equation (4-4.21) as

Diqlm) = f" I w, z)Aq(r, 0dr dt + azfjr [q(r, t) — qpn-(r, z)1Aq(r, 0dr dz
0 r4 0 4
= (w + aT(q - qpri)’ Aq)E = (VJ, Aq)E(4-4.22)
Finally, the equation is in the desired from to identify the gradient of J, which is obtained

from equation (4—4. 1)

84
VJ(r, t, q) = w(r, t) + 0tT[q(r, t) —qpn-(r, t)], (r) on I“4 (44.23)
With the derivation of the gradient VJ complete, the gradient search methods are dis-

cussed.

4-6.0 Optimization Method

The previous equations have presented the derivation for computing the gradient of
function J (q), which is minimized to determine the heat flux. Gradient search methods
are applied to compute the unknown heat flux.

Methods of steepest descent and conjugate gradient are common approaches to com-
pute the search direction. Steepest descent is quite inefficient and not recommended, but is

simple and illustrative. Gradient search methods use iterative techniques to search for the
minimum of an objective function. Beginning with an initial guess q0 for the heat flux, at
subsequent iterations the heat flux is
q"+1(r, t) = q"(r, t) + Aq(r, t), (r) on F4 (4-5.1)
The correction to the heat flux is
Aq(r, t) = —p"p"(r, t), (r) on F4 (4-5.2)
where p" is the search direction and p” is a scalar step size. For the steepest descent

method the search direction is

P"(r, t) = VJ(r, t, q") (4-5.3)
(It is actually the negative of the gradient, but the negative sign is incorporated in equation

(4-5.2).) The conjugate gradient method, which typically converges more quickly, com-

putes the search direction as follows (Arora, 1989)

p" = VJ(r, t, q") + [3"p"‘ 1(r, t), (n at O) (4-5.3a)

85

P0 = VJ(r, t, (10). n = 0 (45.413)

where,

fi’jr {VJ<r. z, ammo, r, q"- l) — We, z. q")1}dr dz

 

B f’j [VJ(rt n-l)]2drdz
t l“4 ’ ’q

___ <VJ(q"), VJ(q"-1)- VJ(q"»E
“won-1)":

 

(4-5.5)

Procedures to calculated the gradient V] (r, t, q") have been outlined previously. Compu-

tation of the step size p" (a positive scalar) remains. The optimum step size is selected
based on reducing J (q) the greatest amount in the current search direction. Equivalently,
the following function is minimized

(MP) = J(q" - Np") (4-5.6)

for (p" 2 O) . The optimal step size, denoted {5", satisfies

(ii = _‘Luqn _ pnp") = O (4-5.7)
do" _ ., dp" .
P" - P P" = P"
The sensitivity problem, which calculates the changes in temperature due to a variation of

magnitude Aq in the heat flux, can be adapted to calculate the change in temperature due
to a variation in heat flux of magnitude p". Then the calculated temperature, which is lin-

ear with respect to q" , can be written as

T(r, z, (1" — pnp") = T(r, z, q") - p"é(r, z) (458)

where 6 is the adapted sensitivity function computed from the problem

V - (k(r)vé(r,z)) = pC(r)[-aa—té(r,t)+v-Vé(r,t)]. (52:11:?) (4-5.9)
o - f

86

-k k—aTB 6,(r z)+h 6(r, z)- _ O , (r) on I}, (1:1’2’3) (4-5.10a)
an (t0<tStf)
_k(r)—a—,é(r, z) = p"(r, z), (r) on F4 (4-5.10b)
3" (t0<tStf)
6(r, O): O, (r) in 52 (4-5.10c)

The adapted sensitivity problem differs from the original sensitivity problem in equation

(4-3.l3) and (4-3.14), only by the right hand side of equation (4-5.10b). In the original
sensitivity problem this is the function Aq. Using equation (4-2.4), the function ¢1(p") in
equation (4-5.6) becomes

<1>(p")- — J(q —p p"): 22f1T<d,,z;"(q —p p ",)>—Y(z)1 dz

j=l

1 ff ,2 _ _ n n 2 2
+20tT zOJnlq (r, t) qpn.(r,t) p p (r, t)] dr dz(4 5.11)

Substituting the variation in the calculated temperature from equation (4-5.8) gives

<1><p")- - 22 f {T(d,,z ;q )-p "91d, z) Y,(z)1 dz

j=l

l I 2
+ 2aT£0Jr,[qn(r, t) — (Jpn-(r, t) — p"p"(r. 1)] dr d! (4512)

After expanding the terms, equation (4-5. 1 2) is rewritten as

up"): 2p 2f161d,.z11dz- p 2ft1m,t;q"1—Y,<z>1é<d,-,z>dz

1:1 j—l
+2 +12]: [T(dj’tzqn )-Y (1)] 2(111+
j=l
l f n n n n n ’1 2
iaTfioIr4{[q (r,t)-CI,,,,-(r, 01-29 p (f,t)[q (r’t)—qpri(rzt)]+lp p 031)] }dr d!

(46.13)

87

Taking the derivative with respect to the step size and solving for the optimal step size [3" ,

as shown in equation (4-5.7), gives

J
2 £:[T(dj, t;qn) — Yj(t)]é(dj, t)dt+ arfiflnpnu, t)[qn(§, t) ‘quz'("z t)]dr dr
—1

 

(Suzi-

J
‘ 2
2 £36”?t)]2dt+arj‘::jr4[p"(r,t)] dr dt

j= 1
(4-5.l4)

Utilizing the relationship between the directional derivative and gradient in equation (4-

3.3) and the directional derivative in equation (4-3.7), the numerator of this function is

J
2 £3“sz t;q") ‘ Yz(’)]é(djz ”‘1’ + arJZjBPWr, t)[q"(~£, t) — Chm-(r, z)]dr dz =
j= 1

(V1 (61"), p">E(4-5.15>

The step size in equation (4-5. 14) is rewritten as
IL. VJ(r, t, q")p"(r, t)dr dt

A ’0 4

p" : :

J J
2 fléwj, t)]2dt + aTflfjr [pn(r, t)]2dr dt 2 "6911., ”"3; aTllpn”:
i=1 0 " “ j=1

(VJ(q"), p">z:

 

 

(4-5.16)
4-7.0 Implementation of the Gradient Method

4-7.1 Whole Domain

The derivation to solve the IHCP using a gradient search method and adjoint equation
approach have been presented. A summary of the solution with these methods, in the form
of an algorithm for a whole domain implementation are given next. Equations presented in

this section are identical to those derived earlier and are summarized here for clarity.

88

1) Select initial values for the function q°(r,t) (n = 0). A common choice is
q0(r, t) = 0. If prior information is available set an-(r, t).

2 i) Solve the direct problem for the temperature T(r, t;q")

V-(k(r)VT(r,t)) = pC(r)[—a—T(r,t)+v-VT(r,t)], (r) i" 52 (4-6.1)
at (10 < t s. If)
_k, inr, z)+h, T(r, z) = f2.(r, z), (r) on 13"“ =1’2’3) (4—6.2a)
3" , (tO <15 1f)
-k(r)3.-T(r. t) = q"(r, t) . (r) on 1‘4 (4-6.2b)
3" (r0 < z g If)
T(r, 0) = T0(r) , (r) in Q (4-6.2c)
2 ii) Evaluate the residuals
eJ-(ti) = Yj(t,-) — T(dj, ti;q") (4-6.3)
3 i) Solve the adjoint problem to determine w(r, I)
a J
V . (k(r)V\1I(r, z)) + [pC(r)§—t-w(r, z) + v - V[pC(r)w(r, z))] + X ej(t)5(r—-dj) = 0
j: 1
(r) in Q (4-6.4)
(r0 < t S If)
-k.- 587w. z) + 1h.-— pC(r11v-zz)1w(r.z) = o, ”’ °" F" (' z" 2’ 3) (4-65a)
n (10 < z g If)
k(r) give, z) + pC(r)<v - zz) w. z) = 0. m °" r4 (4-6.51)
" (tO <13 1f)
w(r, t)|t = z, = O, (r) in Q (4-6.5c)

3 ii) Evaluate the gradient VJ (q")

89
VJ(r, t, q") = w(r, z) + ozT[q"(r, t) —qpn-(r, t)], (r) on F4 and (t0< t 3 if) (4-6.6)

3 iii) Compute search direction p"

p"(r, t) = VJ (r, t, q") - Steepest Descent Method (4-6.7)
p"(r, t) = VJ(r, t, q") + [3"p” ‘1(r, t), (n ¢ 0) Conjugate Gradients (4-6.8)
p0 = VJ(q0), n = O Conjugate Gradients (4-6.9)

where,

fifjr {VJ(r’ I, qn)[VJ(r’ t’ q") - VJ(I’, t, (In— l)]}dr dt
fin = 0 4

 

2":er [VJ(r, t, q"'1)]2dr dt
___ (WW), VJ(q"")—VJ(q")>zs

 

 

 

2 (4-6 10)
IIVJ(q"-1)Ilz:
4 i) Solve the (adapted) sensitivity problem for 8(5, t)
V - (k(r)Vé(r, t)) = pC(5)[-a—é(r, t) + v ~ V60, 1)], (r) 1n Q (4-6.l 1)
at (10 < t _<_ tf)
(i =1, 2, 3)
—k2. 23—2260, 1) + h,- é(r, t) = O, (r) on I",- (4-6.12a)
(r0<ts If)
_ F
~k(r)—Q;9(r, t) = p"(r,t), (r) on 4 (4-6.12b)
3" (to <13 If)
é(r,0) = 0, (r) in Q (4-6.12c)
4 ii) Compute the optimal step size (1
f
‘7 n n
An _ 10L} J(r, t, q )p (r, t)dr dt — (VJ(q"), Pn>E
p - J 2 2 _ J 2 2
f - I n - n
.lejulewj, 1)] dt+a7£0jr4lp (r, 1)] dr dt _21I|9(d,-.t)llg+azllp HE
J : j :

(4-6.l3)

90
5) Compute the improved value for the heat flux
q" + '(r, t) = q"(r, t) — fi"p"(r. t). (r) on F4 (MM)

6) Check for convergence of the estimated heat flux

2
ff],- [qnflvz t)-q"(r, 1)] 4’ dt = Ilq”+‘-q"||:~<€ (4-6.15)
‘0 4

If convergence has not been obtained let n = n + 1 and return to step 2 with the updated
heat flux. If convergence has been met, use the residual principle to verify that the correct

magnitude of (IT was used. The residual principle is

J
15(q") = 2 J:’[T(d,-, q")—Y,-(t)]2a't=1:52 (4-6.16)
1:1 °

2 . . . . .
where 5 IS the expected no1se or error 1n the measurements and 1‘ 2 1 1s a relaxat1on
parameter.

An iterative regularization method can be obtained in the above algorithm by allowing

the regularization parameter (IT to go zero. However, convergence is not evaluated by

checking that by changes in the heat flux are small as in equation (4-6.15). Instead, in iter-
ative regularization, the residual principal is used to select the number of iterations such
that equation (4-6.l6) is satisfied and the sum-of-squares function is reduced to its
expected value; once this is achieved the iteration process is stopped.

Sequentially implementing this method is discussed next. Insight to the solution pro-
cedures using the gradient method with adjoint equation approach is given in Section (4-

8.0) at the end of this chapter.

9 1
4-7 .2 Sequential Implementation

The usual application of a gradient method with adjoint equation approach would
solve the three problems over the whole time domain (tO < t S t f) , obtaining the estimated
function q(r, t) for all time. Consequently, any inaccuracies in the estimated function at
later times might influences the estimated function at early times. This dependence results
from simultaneously estimating the function for all time.

An alternative approach uses the sequential idea applied in function specification
method. A sequential approach is quite effective for estimating the heat flux. In addition, it
possesses other positive attributes. The components that are estimated at later times are not
influenced by the early estimates, less computer storage (and possibly time) is required,
and it is possible to linearize over the future interval and more efficiently solve non-linear
problems. Benefits of a sequential method are discussed in more detail in Section 5-2.0.

The sequential procedure can also be implemented for a gradient method using an
adjoint equation approach. This procedure applies the method over smaller time intervals.
Assuming that the function q(r, t) is known for times t S rm _ 1 , the unknown heat flux
over r-future time steps, tm _ 1 < t S rm + r_1 is estimated. Equations for the whole domain
implementation are identical to those for the sequential implementation; however, they are
solved over the time interval tm _ 1 < t S tm + r _ 1 . After obtaining a solution, the time inter-
val is shifted by one (or possibly more) time step(s) and the procedure is repeated until the
whole time domain is covered. In contrast to function specification it is not required to
assume a functional form for the heat flux, but such a procedure may be beneficial. With
this sequential approach, many shorter duration problems are solved, compared to the

whole domain approach that solves one large time duration problem. Operationally. it is

92

straight forward to use gradient search methods in a sequential manner. Implementing
such a method, and the associated practical issues, are not as straight forward. Chapters 5
and 6 study and provide insight to this implementation.

It is anticipated that this sequential approach will result in great reductions in the use
of computer memory and have the potential to substantially reduce computational time.
More savings are expected as the number of measurements become very large. These
potential benefits are to be explored and are important parts of the proposed research. A
more in-depth discussion of the benefits associated with the sequential approach is given
in the subsequent chapters (Section 5-2.0).

The formulation of the problem and the describing equations are nearly identical for a

sequential implementation. However, the sequential implementation solves a series of

smaller problems over the time domain. The time domain (to < t S t f) is broken into a
series of sequential intervals each of length (tm_l < t S tm + r_ l) and solved for

m = l, 2, . . . .M, where M = (tf- t0)/At— r, AI is the time step, and r is the number

of future time steps. Assume that the heat flux for the sequential intervals up to and includ-

ing time tm _ 1 is known. It is desired to compute the heat flux for subsequent times. Simi-

lar to the whole domain implementation, the solution is presented in the form of an

algorithm.

1 For the first sequential interval (m = 1) and iteration (n = 0) select the initial values for
the function qO. Future sequential intervals (m > 1) use the estimates at the previous

sequential interval to specify the initial value. A common choice is q0(r, t) = 0. If prior

information is available set qpn-(r, t) .

93

2 i) Solve the direct problem for T(r, t;q")

V - (k(r)VT(r, z)) = pC(r)[3T(r, z) +v . VT(r, n], (r) i“ Q (4-6.l7)
at (tm_1<tStm+,_1
_k, ainr, z)+h, T(r, z) = f2.(r, z) , (r) 0" Ff“ =1’2’3) (4-6.18a)
n (tm—l<tStm+r-—l)
—k(r)-:TT(r, z) = q"(r, z) , (r) 0" 1‘4 (4-6.18b)
n (tm—l<tStm+r—l)
T(r, tm_]) = T(r, tm_1), (r) in Q (4-6.18c)

2 ii) Evaluate the residuals

eJ-(t) = Yj(z)—T(dj,z;q") ,(tm_,<tszm+,_l) (4-6.l9)
3 i) Solve the adjoint problem to determine w(r, t)
a J
V . (k(r)Vw(r, r))+[pC(r)§7w(r,t)+v - V[pC(r)w(r, z)]]+ 2 ej(t)5(r—dj) = 0
j=1
(r) in Q (4620)
(tm-l <ts‘m+r—l
F, '=l,2,3
-k.-ierthz-parxv'fin Wat) = 0. (”on ‘ (I ) <4-6.21a>
an (tin—l<tgtm+r—l)
F
k(r) inl(r,t)+pC(r)(v-fi) Wm) = 0. (r) on 4 (4-6.21b)
an (rm—l<t'<‘tm+r-l)
w(r, 0|! :2 = 0, (r) in Q (4-6.2|c)
3 ii) Evaluate the gradient VJ ( q")
F
VJ(r,t,q”) _—. 111(r, t)+0tT[q"(r, t)—q2m-(r, 1)], (r) on 4 (4-622)
(tin—l <tstm+r—l)

94

3 iii) Compute search direction p"

p"(,r t): VJ (r t, q")- Steepest Descent Method (4-6.23)
P"("z t) = VJ(r, 1. q") + [3”p" ‘1(r, t), (n ¢ 0) Conjugate Gradients (4-6.24)

p0 = VJ (qo), n = O Conjugate Gradients (4-6-25)

where,

[::” “'n[ {VJ(r, z, q")[VJ(r, z, q")— VJ(r, z, q" 1)1}dz- dz
B" =

 

[:m _ [r [VJ(r, t, q"‘1)]dr dt
= (VJ(q"), VJ(q"' 1) — VJ(q")>E

 

 

 

2 (4—6.26)
IIVJ(q"“)||E
4 i) Solve the adapted sensitivity problem for 6(r, t)
V (k(r)V6(r, t))= pC(r)|:(%- —6(r, t)+v V9(r, 0], (r) 1n 0 (4-6.27)
(rm—l <tStm+r—l
- - F, ' =1, 2, 3
—k 1,90,04,19. 9(r,t) = 0, (r) on ' (l ) (4-6.28a)
an (tm—l<tStm+r—l)
(r) on F4
-k(r)§- "6(r, t)- — p"(r, I), (4-6.28b)
(tm—l<tStm+r—l)
6(r, 0) = O, (r) in Q (4-6.28c)
4 ii) Evaluate the optimal step size 6
[JIM IJF VJ(r, t)p"(r, t)dr dt
A 1 —1 4
P": m
it: '[é(d 1)]2dt+01 [”‘”'[ 1 "(r t)]2dr dz
I j, T (In-l 1.‘4 p ,
j-
V] n 2 n ‘
= J < (61 ) P >15 (“529)

z ||é<d,. 0112+ odnpnnz

j=l

95

5 Compute the improved value for q

F
WHOM) = CI"(r,t)-15"p"(r,t). (r) on 4 (4-6.30)
(tm-l<tStm+r-l)

6 Check for convergence of the estimated heat flux

J'r JJm+I-l[(qn+l(r, t)_q"(r, t))]2dt dr ___ llqn+l_qn”2<8 (4-6,3])
4 ’m—1

If convergence has not been obtained let n = n + 1 and return to step 2 with updated heat
flux. If convergence has been met, store the heat flux for sequential analysis interval (store
l-components)

qes,(r, z) = q" +1(r, t), (zm_ 1 < IS tm+1_1) (4-6.32)
Then index sequential interval m = m + I and reset iteration to n = 0. Use estimated
heat flux from the current sequential interval to partially set initial guess for next sequen-
tial interval (an additional I-components must be specified to complete the initial guess).

Return to solution of direct problem (step 2.)

7 After final sequential interval check that proper magnitude of (1T was used with the

residual principle
J 2
z 2
J(q") = 2r [T(dj, q”) — Yj(t)] dt =- 15 (4-6.33)
!
1= 1 °
2 . . . . .
where 5 1s the expected no1se or error 1n the problem and t 2 1 1s a relaxation parameter.
The selection of the Tikhonov parameter for this case does not take advantage of the
sequential implementation to improve the selection of the parameter during a solution.
Selecting the parameter after completing the solution for the whole time domain was cho-

sen to allow a direct comparison between a sequential and whole domain solution. It is

96

possible to sequentially select the magnitude of the Tikhonov parameter at each sequential
interval. In a sequential selection of the parameter, equation (4-6.33) is evaluated at each
sequential interval and criteria are specified to vary the magnitude of the Tikhonov param-
eter. In a sequential application of equation (4-6.33) the results would likely be more sen-
sitivity to local errors, and some constraints based on previous values of the Tikhonov may
be need to prevent large fluctuation in the parameter from one sequential interval to the

next.

4-8.0 Parameterizing the Solution (Finite-Dimensional Problem)

The derivation up to this point does not assume information concerning the functional
form of the unknown heat flux. The problem, in general, has been considered infinite
dimensional. This approach permits estimating a functional representation of the heat flux
without having to specify information concerning the unknown function. Of course for
realistic problems, obtaining a solution requires a numerical procedure. Any numerical
procedure chosen will require discretizing, both on space and time, and consequently
making assumptions about the form of unknown function. These assumptions are not too
restrictive for most cases and the function can be estimated at the discrete grid points. If
however, information is available, or assumed, concerning the heat flux, the solution can
be constructed as a special case of the generalized solution methodology presented for the
infinite dimensional problem.

Consider the case that information about the form of the heat flux is known. A general
form for the known information can be represented as

P M

q(r,!) = 2 zqz,k¢z(t)<pk(r) (4-71)

k=1i=l

97

where ¢z(’) and (pk(r) are basis functions on time and space and qi’ k is the parameter-

ized heat flux. When the heat flux is represented in the form of equation (4-7.1), solving
for an adjoint function is not required to compute the gradient VJ . Standard differential
calculus can be used to compute the partial derivative with respect to each component of

heat flux BJ/qu, k. When there are a relatively small number of components this approach

is suggested. However, for the multi-dimensional problem a large number of components
are required to represent the spatial, time varying heat flux. Hence, the solution for the
finite dimensional representation of the heat flux is presented as a special case of the infi-
nite dimensional problem.

As stated previously, solving realistic problems requires using a numerical procedure.
Regardless of the numerical method selected (finite element, finite difference, boundary
element, or finite control volume) the spatial and temporal domains required discretizing.
Consequently, even the infinite dimensional problem requires the heat flux to be
represented in a finite dimensional form, as shown in equation (4-7.1), when solving with
a numerical procedure. For the finite control volume numerical method, assuming that the

quantities are constant over the control volume, P is the number of control volumes on

surface F4 and M is the number of time steps. It is well understood that improved

numerical accuracy is obtained by decreasing the size of the control volumes and time step
when solving the direct problem. Although this improves the numerical accuracy of the
direct solution, it also increases the number of spatial heat flux components and decreases
the magnitude of the time step. Estimating a larger number of spatial components or at

smaller time increments makes the inverse problem more ill-conditioned. Hence, refining

98

the spatial and temporal grid size for a numerical method will not necessarily improve the
accuracy of the inverse problem, in fact it will likely produce a larger variability in the
estimated function. If however, some smoothing or prior information is used, as in
equation (4-7.1), accuracy will not be significantly affected.

Modifications of the solution for a heat flux represented as in equation (4-7.l) are
mainly in the computation of the gradient. Implementing the finite dimensional case using

an adjoint equation approach is given in Jamy et al. (1991) and Alifanov (1994). Assum-
ing the update to the heat flux, Aq(r, t) has the same form as the heat flux, and substitut—

ing these approximations into the directional derivative in equation (4-4.22) gives

P M
DA,J(q) = 2‘, Z J” J we, z)¢,-(z>zp,.<r>Aq(r, z)dr dz
1' I F,

k = 1:
P M P M
z
+ az 2 2 2 2 J0 Jr41¢.(z)<p,.(r>q(r. z) — q,.,-<r. z)1¢,(z)<p,(r)Aq(z-, t)dr dz
l=lj=lk=11=1
(4-7.2)
Equation (4-7.2) can be arranged to determine the gradient
P M
DAqJUI) = Z 2 VJ z, kqu, J. (47.3)
k=li=1

where

P M

V‘IiJz = Yz,k+aTZ 2Si,k,j,1[q(rzt)—qprz(rzt)]~ (1:1’2’ W’M) (4-7.4)

I=Ij=l (k:1,2,...,P)

P M
7.1. = 2 2 f ’ J we. z)¢z,-<z><p,.(r)dr dz (4-7.5)

k=li=l “F.

Si,k,j,l = J:In¢,-(t)(Pk(r)¢,-(t)(p,(r)dr dz (4-7.6)

99

Thus, once the adjoint function \JJ(r, t) is available, 72-, k is computed by integrating with

the known basis functions. Similarly the weighting for the regularization S i, k, j. I is com-
puted by integrating.

Restricting the more general infinite dimensional problem to a finite dimensional
domain is straight forward. After solving for the adjoint function, the gradient for the finite
dimensional representation of the heat flux is an integration of the adjoint function with
the associated basis function. Other quantities computed in the solution, such as parame-

ters for the search direction and step size, are performed in the finite dimensional space,

 

i.e.,( , )R, and II IR”.

4-9.0 Summary

A formulation to solve the multi-dimensional IHCP using a conjugate gradient search
method with an adjoint equation approach was derived in Section 4-1.0 to 4—6.0. The
describing differential equations for solving in a whole domain implementation, herein
referred to as the whole domain gradient method (GM), were presented in algorithm for—
mat in Section 4-7.l. The formulation for a method that solves the IHCP sequentially in
time was also given in an algorithm format in Section 4-7.2. Studying this sequential
implementation, herein called the sequential gradient method (SGM), is a major objective
of this dissertation. The algorithms, for both GM and SGM, were formulated for the most
general case, allowing the estimated heat flux to be described as a function (infinite dimen-
sional). As a special case, restricting the estimated heat flux to a finite dimensional space

was discussed in Section 4-8.0.

100

Solving the IHCP with a gradient method requires solving three problems, which are
denoted the direct, adjoint, and sensitivity problems. Results from the solution of each
problem provides an input for the subsequent problem. Residuals computed from the solu-
tion of the direct problem, equations (4-6.17) and (4618) are the driving term in the
adjoint problem in equation (4-6.20) and (4-6.21). The search direction computed with
information from the adjoint function is the driving term in the sensitivity problem in
equation (4-6.27) and (4-6.28). If the measured temperatures are known “exactly” i.e.,
residuals are zero, all three problems have a homogeneous solution. For inexact data, as
iterations progress the solution of the three problems should converge towards a homoge-
neous solution.

The remaining chapters study the solution of the IHCP using the GM and SGM. In
Chapter 5 a detailed discussion of the benefits of a sequential implementation are
addressed. Because many of these benefits are related to the numerical solution, a discus-
sion of the numerical aspects are given first. In the remainder of Chapter 5 the methods are
applied to the one-dimensional IHCP. Although a two-dimensional formulation is given, a
one-dimensional solution is obtained as a special case of the two-dimensional formulation.
In general, it is not recommended that this solution technique be applied to a one-dimen-
sional problem. However, to gain insight to the SGM and understand it better, a one-
dimensional solution is an appropriate beginning to the investigation. After investigating
the one-dimensional problem, Chapter 6 progresses to study the two-dimensional prob-

lem.

Chapter 5

APPLICATION OF THE SEQUENTIAL GRADIENT METHOD:
ONE-DIMENSIONAL IHCP

5-l.0 Introduction

Cases are presented in this chapter to study the sequential gradient method. The accu-
racy of the sequential implementation, as compared to the whole domain, is also quanti-
fied. Addressing issues of accuracy for inverse methods is an important topic. Though
there are many methods for solving the IHCP, few methods have had detailed studies con-
ducted on their accuracy. Investigations of the accuracy of inverse methods are presented
by Beck and Murio (1986), Raynaud and Beck (1988), and Scott and Beck (1989). Beck
and Murio (1986) compared a whole domain solution with a sequential implementation
combining function specification and Tikhonov regularization. The work by Scott and
Beck focussed on the sequential regularization method while Raynaud and Beck com-
pared four methods for solving the IHCP. In these papers, however, simulated cases, with
known heat flux, are analyzed to quantify the accuracy of the inverse methods. These ideas
are also applied in this investigation. In particular, standard test cases for which the heat
flux is known are analyzed and the accuracy is quantified for the sequential implementa-
tion as well as the whole domain approach, which is the usual application of this method.
Hence in addition to studying the sequential method, this dissertation provides insight to

the whole domain solution.

101

102

An important issue in an inverse analysis is the effect of measurement errors on the
solution. Many researchers study this effect modeling the errors with standard statistical
assumptions, which is beneficial and also pursued in this investigation. Additional insight
is gained from an analysis with actual experimental measurements. Hence, in addition to
the studies with simulated errors, an experimental case is presented. The case has the
unique benefit, in that transient temperature and heat flux histories are both measured.
Since the heat flux is measured, the “true” heat flux, within some experimental uncer-
tainty, is known. This experiment permits comparison between the estimated and mea-
sured flux.

A new competitive sequential method is presented in this chapter. In addition to
studying the sequential gradient method a modification of the sequential implementation
is suggested that improves the computational requirements of the sequential method. The
new method combines function specification with the sequential gradient method. Imple-
menting a gradient solution, but assuming the heat flux is a prescribed function over the
sequential interval is studied. Prescribing a functional form for the heat flux over the
sequential interval is shown to reduce the computational requirements compared to a stan-
dard sequential implementation.

Associated benefits of a sequential implementation of gradient methods and the
numerical procedures to obtain a solution are discussed in next section. Several test cases
to simulate experimental data are given in Section 5-3.0. Results of the one-dimensional
analyses with simulated exact data, data corrupted with measurement errors, and experi—

mental data are presented in Section 5-4.0.

103

5-2.0 Advantages of a Sequential Method and Numerical Solution of the
IHCP

This section discusses the benefits and advantages of a sequential implementation for
solving the IHCP. Because these advantages are closely related to the numerical solution,

the numerical method to solve the IHCP is discussed first.

S-2.1 Numerical Solution

A numerical solution for the two-dimensional IHCP is developed using a finite con-
trol volume (FCV) method. Numerically solving the IHCP using a gradient method and
adjoint equation approach requires the solution of three partial differential equations rep-
resenting the direct, adjoint, and sensitivity problems. Describing equations formulated for
the three problems are given in equation (4-6.17) and (4618), equation (4-6.20) and (4-
6.21), and equation (4-6.27) and (4-6.28), respectively. Fortunately, the differential equa—
tions are similar and numerical computations required for one problem apply for the other
problems. A detailed derivation of the FCV equations are given in Appendix A.

The numerical solution was developed for the two-dimensional problem using an
alternating direction implicit (ADI) scheme. A one-dimensional solution is obtained as a
special case. The FCV equations for the direct problem are represented as a two step solu-
tion

1K11{T""”2} = {0"} (5-2.1)

} = {EMW} (5-2.2)

n+1

[K 21{T
where [K 1] and [K2] are the standard FCV tri-diagonal matrices, {T"} is a vector of

n+l/

. 2
unknown temperatures for t1me n, and {0"} and {E } are vectors of known

104

information for the respective time steps. The two sets of equations represent the two steps
of an ADI scheme.
All three problems, direct, adjoint, and sensitivity, require solving a set of equations

similar to those shown in equation (5-2.1) and (5-2.2). Since the problems are similar, the

tri-diagonal matrices ([K1] and [K 2]) are identical for the three problems and require

computation only one time for constant thermal properties, i.e., a linear problem. Further-
more, when solving the tri-diagonal set of equations, computational savings are possible
for the linear problem. This is true for a whole domain as well as sequential solution.
However, for non-linear problems with a sequential solution it is possible to temporarily
linearize and benefit from the computational savings associated with linear problems.

Aspects of this procedure are discussed below.

5-2.2 Sequential Solution of IHCP

Benefits of applying a sequential solution for the IHCP are addressed. Several sup-
porting reasons for a sequential solution are discussed. Although not all reasons are inves-
tigated in this dissertation, nor verified, it suggests the possibility of growth for the

sequential gradient method.

1. A sequential method can be implemented in an on-line or “real” time mode. For exam-
ple, in the monitoring of surface heating during the flight of a space vehicle, data can be

collected for a short period then the surface heat flux can be computed in “real” time.

2. Estimated values much later in time do not affect the estimates at early times when
solving sequentially. When solving in a whole domain approach, all temporal compo-
nents are simultaneously estimated. Consequently, errors in estimated components will
pervade all estimates. It is physically unreasonable that this occurs because of the para-

bolic nature of the describing equations. Furthermore, values near the final time, which

105

are known to be difficult to estimated, must be estimated and the errors in these compo-
nents may impact the others. Since a sequential approach solves on a subset of the
whole time region, or sequential region, estimates at later times do not impact early

estimates, which is physically reasonable.

. For nonlinear problems, due to temperature dependent thermal properties or other tem-
perature dependent coefficients in the differential equation, the sequential method per-
mits the problem to be temporarily linearized. This is an important point, because it is
not possible when solving on the whole time domain. Since the problem is solved on
the shorter sequential interval, it is valid to assume that the temperature dependent
quantity is constant during the sequential interval (r-future time steps). The validity of
this assumption depends on the degree of non-linearity in the problem. In most applica-
tions the non-linear term is not known with sufficient accuracy to justify accounting for
the non-linearity during the temperature change seen in a sequential interval. To linear-
ize, thermal properties and temperature dependent variables are evaluated at the initial

temperature distribution of the sequential interval. Consequently, computing matrices

[K1] and [K 2] is required only at the beginning of each sequential interval. The matri-

ces apply for all three problems and it is not required to recompute the matrices during
the iteration procedure. Conversely, the whole domain approach requires re-computing
the matrices at each time step as the temperature varies over the time domain. Matrices

must be recomputed at every iteration.

. Related to point 3., is the computational savings that can be realized in the solution of
the tri-diagonal system of equations. Since the problem is linear (temporarily) for a
sequential approach, the intermediate coefficients required to solve the tri-diagonal sys-
tem can be saved after solving the direct problem and need not be computed until the

next sequential interval, when the initial temperature, and coefficients, change. Examin-

ing the algorithm illustrates the point. Assuming the tri-diagonal matrix [K1] has main

diagonal C1,, sub-diagonal A1,,- and sup—diagonal B the system is solving by the

1,1"

Thomas algorithm using the recursive relations (Anderson et al., 1984)

B0 = C1, 0 (5-2.3a)

 

A1 iBl z—1 .
51: 12— ' ' ,z=1,2,...N (5-2.313)
13,--1
n
_ Do
yo — — (5-2.3c)
Bo
D'.'-A . .
y,- = ' "‘Y“l,i=1,2,...N (5-2.3d)
B.-
The solution of the temperature is
1/2
T2: = 7N (5-2.4a)
B _Tzz+1/2
Tj””2=y,——"—'—B”—‘,z=N—1,N—2,...1,o (5-241))
1'

Since the problem is linear (temporarily), the coefficients A 1’ ,- , C1, 2. , and, 3U are the

same for the direct, adjoint, and sensitivity problems. Consequently, after computing

and storing [3,. in equation (5-2.3b), this step is by-passed in subsequent solutions. A

similar procedure is used to store the coefficients from the second step in the ADI
scheme. Computations in the following steps, equation (5-2.3c) and (5-2.4a) and (5-
2.4b), are required for every solution. For a linear problem, both the sequential and

whole domain methods can benefit from this computational saving.

. Solution of the adjoint and sensitivity problems can use the zero initial condition and
fact that driving terms in the equations tend to have oscillating signs to reduce computa-

tions by marching on space and/or time. The possibly oscillating signed driving terms,
residuals e 2(t) in the adjoint problem and search direction p"(r, t) in the sensitivity

problem, will localize the effects. It may be possible to consider the solution on a subset
of the time domain or physical domain because the effects are localized. A whole
domain approach might benefit from such an approach, but because the sequential

implementation focuses on a restricted time region it is more likely to benefit.

107

6. The magnitude of the Tikhonov regularization parameter can be adjusted in a sequential
implementation. Again, this is not restricted to the sequential implementation, but is
more practical to implement in a sequential method because current information is
available to select the magnitude of the Tikhonov parameter. This permits more accu-

rately estimating rapidly changing functions with less bias.

7. Additional heat flux components can be retained in a sequential implementation. This is
an advantage compared to a standard function specification solution, because additional
bias is not introduced with the function specification approximation. More components

from the solution over a sequential interval therefore may be retained.

5-3.0 Simulated Test Cases

Several standard one-dimensional cases are examined. Test cases include heat flux

functions that vary in a triangular, step, and sinusoidal manner with time. The one-dimen-
sional geometry is a slab of thickness L with an unknown heat flux at x = O , insulated on
the surface at x = L, uniform initial temperature, and constant thermal properties. Since
the IHCP is very sensitive to measurement errors, to reduce the possibility of numerical

errors, analytical solutions (when possible) are employed to generate simulated data.

Mathematically the problem is represented as

 

 

 

a Tj.’ aTj'
2 = ‘—+ (5-3.I)
ax“ az
ar‘.‘
_ a i = q:(z*) = ? (5-3.2a)
x 22+:
6T:
2 = 0 (5-3.2b)
ax

T:(x+, 0) = 0 (5—3.2c)

108

where i = 1, 2, 3 for the three prescribed heat fluxes studied. To remove the dimensional
dependence and facilitate the investigation, dimensionless variables were introduced in

equation (5-3. 1) and (5-3.2). The dimensionless variables are

 

q: = 3- (5-3.3)
4N
T--T
TI = ' 0 53.4
' (qNL/k) ( ’
1+ = @1291 (53.5)
L
+ x
= — - .6
x L (5 3 )

Symbol q N is a nominal value of the heat flux, k and pC are the thermal properties, T0 is

the uniform initial temperature, and L is the plate thickness in equation (5—3.3) through

(5-3.6).
Three forms of the prescribed (assumed unknown) heat flux (q?) are examined. A tri—

angular heat flux

0 (I+<0)

+ < +< +
CIR?) = [ t (0” _Ih/z) (5-3.7)
t+—2(t+-—I;/2) (tZ/2<I+SI;)

O (t+>t;)

 

K

which increases linearly to a magnitude of t;/2 then decreases linearly until reaching

zero, and is zero elsewhere. The second variation is a step change in the heat flux

109

 

0 (t+<0)
+ +
‘12“ ) = i 1 (OSt+St;) (5-3.8)
J O (t+>t;)

which increases to a magnitude of one at t+ = 0 , remains constant until t+ = t; and then

returns to zero. The final heat flux is a sinusoidal heat flux. It represents one-half the

. 1t 31: . .
perrod {—5, 3-) of a s1ne wave between zero and t+ = I; , and 1s zero elsewhere

 

,
O
(t+<0)
+z+ — 1 - ’+ 1 + + 539
‘13()-< isrn—n2:+§ +1 ogzgzh ('-)
’zz
[ O (t+>t;)

A constant is added to the sinusoidal heat flux to maintain a nonzero flux.

Analytical solutions are readily generated for these three cases. Solutions can be
obtained using Greens functions (Beck et al., 1992). Although the test cases mentioned are
standard cases and it is not difficult to solve for the temperature, closed form expressions
are not typically given in the literature. Beck et al. (1985) give the closed form solution for
a linear heat flux at the boundaries (x = 0, L). For completeness, the closed-form expres-
sions for the temperature solutions with the three specified heat flux functions are pro-
vided.

For the triangular heat flux the solution is built using superposition of a heat flux

which increases linearly with time (q+ = t+ ). The temperature due to the linearly varying

heat flux is

+2

+++_t 1+2+l 1M1+31+21
¢l(x,t)——é—+[x —X '1' “It +[2—4x ~6x +6x —'4-5]

+2 0° 1 2 2+
+7212 —4c cos(mrtx )[exp(——m 1: t )](5-3.10)

The temperature solution for the triangular heat flux in equation (5-1.7) follows as

 

u*<0)
+ + + + +
.x , t <: 3<
41m". t“) — 2<1>I'(zc",t+ — III/2) (z;/2 < z+ _<. z;)
1 (q(x", 1+)— 2¢f(x*, z+ — z,j/2) + ¢f(x*, z” — z;) (,+ > ti)

(5—3.l 1)

Similarly, the temperature solution for the step help flux is built using superposition

of a constant heat flux.The temperature solution for a constant heat flux (q+ = l ) is

00

¢;(x+, t+)= t + [1x +2 —x+ + l] — —2— —1—2- cos(m1rx )exp(—m211:2 t ) (53.12)
2x 3 K2 17712
n =

For the step change in the heat flux, equation (5-3.8), the temperature solution is

(1+ <0)
T+ + + — 4 + + + + +

 

+ + + + + + +
\ ¢2(X 3t )—¢2(X rt _th) (t+>t;)
The sinusoidal heat flux in equation (5-3.9) is the sum of a sine function and a constant,

hence, superposition of these solutions is used. The sinusoidal part of the heat flux,

q+ = sin [—1t(21+/t; + l/2)] , has a temperature solution

111

+ 00
th t+ l 1
¢;(x+, f) = ECOS[—n[2-+ + 2]] +4t; 2 [m4 4 +2 2:|cos(m11‘.x+)

th m=1 “’1. +41:

+ +
{21tcos[—1t[2-t: + a] + mzrtzt;sin [—1c(2t—+ + g] + mZTtZIZexp(—m21t2t+)} (S-3.14)
th ’11

The temperature solution for the sinusoidal heat flux in equation (5-3.9) is

 

 

1 0 (t+ < 0)

73oz“, z“) = i1 130*, z“) + ¢§(x+. z”) (0 S ti 5 ’Z)
1 ¢§<x+, z“) — «1306', z+ _ z;) + ¢§(x*, z*) — <1>§(x+,t+ — 12) (z+ > t2)

(5315)

Data for the simulated temperature and heat flux are shown in Figure 5-1, 5-2, and 5-
3 for the triangular, step, and sinusoidal heat flux. All three cases demonstrate the effect
that makes the IHCP difficult, mainly the lagged or delayed response in the temperature at
the back surface due to the heat flux on the opposite surface. The temperature response is
delayed relative to the application of the surface heat flux. The more gradual the variation
in the surface heat flux is, the greater this delay becomes, as seen by comparing Figure (5-
2) and (5-3). The temperature response at the sensor location is damped also. That is, the
magnitude of the temperature response at the sensor location is smaller that the response at
the surface where the heat flux is estimated. These temperature data are used to estimate

the surface heat flux.

112

0.6 I I F

 

0.5 - _

0.4— q mm

 

I

0.3 A

0.2 r 2

Temperature / Heat flux (dimensionless)

 

 

 

Figure 5-1 Simulated temperature and heat flux data for the one-dimensional triangular
heat flux

113

1.2 I I I

 

 

 

0.4 -

I

0.2

Temperature / Heat flux (dimensionless)

 

 

 

 

 

 

1.5 2

Figure 5-2 Simulated temperature and heat flux data for the one-dimensional step heat
flux test case

114

 

 

 

 

 

 

 

1 I I T I
A 0.8»-
m
as
.5 T*(1,:*)
<0 _ +
E
3
E?
a: 0.4L
2;;
Q)
m
\
8 0.2“
D
a
H
0)
Q.
E o_
[—
_O.2 1 l l L
-0.5 0 0.5 1 1.5

Figure 5-3 Simulated temperature and heat flux data for the one-dimensional
sinusoidal heat flux test case

115
5-4.0 1D Results - Simulated Measurements

One-dimensional test cases are presented to study the sequential implementation.
Implementing such a method is not recommended for the usual one-dimensional prob-
lems; it is applied in this case to develop an understanding of this new method. (Nonlinear
problems may be an exception for using the sequential method for one-dimensional prob-
lems.) Three functional forms of the heat flux, triangular, step, and sinusoidal, serve as test
cases. Investigations with exact data and data corrupted with errors are conducted. Errors,
denoted 8n(dj, t) , are assumed additive

Y(dj, t) = T(dj, t) + en(dj, t) (54.1)
with a normal distribution, zero mean, constant variance, and uncorrelated. The true tem-

perature T(dj, t) is computed from the analytical solutions, see Section 5-3.0. Random

errors are generated with routines in MATLAB.

5-4.l Exact Data (no Measurement Errors)
The three test cases are analyzed with exact data (an = O) for a dimensionless time

step of 0.06 and a total of thirty time steps. Values selected for r are five, eight, and ten.
Results are shown for the triangular heat flux in Figure 5-4. Graphic results are not shown

for the step and sinusoidal test cases using exact data, but the analyses are discussed. Esti—
mates for the final sequential interval (t f — rAt < t < If) at the end of the time region are
not considered reliable and not reported. The results of the estimations are similar with

less than a five-percent difference in the whole domain and sequential implementation;

errors are largest near areas of sudden change in the surface heat flux.

116

 

 

Heat flux (dimensionless)

 

 

 

0.8 ' T T T I I t I . I
o——oGM 4
H SGM r=5
0 6 _ <>"—">SGI\A r=8 ..
' H SGM r=10
— Exact
0.4 - -
0.2 - .
0.0 " i-T "3‘: ;;:: _
-02 . I 1 1 1 1 . 1 1 1 .
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0
+
t

Figure 5-4 Estimated heat flux for one-dimensional triangular test case with exact data

117

Table 5-1 Estimation results for exact simulated data

 

 

 

 

 

 

 

 

Iterations
. Tikhonov Comp S 1’ 0 D
Analysrs per time 2
Method domain 0‘7" total seq int (sec) (°C) W/m

 

 

Triangular Flux (O’ = 00°C)

 

 

 

 

 

GM whole 0.001 7 - 2.14 8.59 13-04 0.0094
SGM r = 5 0.0001 1 78 3.0 3.64 8.64 E-04 0.0046
SGM r = 8 0.00087 87 3.8 6.86 8.72 E-04 0.0094
SGM r = 10 0.00107 87 4.1 8.64 8.72 E-04 0.011

 

 

 

 

 

 

 

 

 

Step Flux(6 = 00°C)

 

 

 

 

 

GM whole 0.001 8 - 2.38 3.12 E-03 0.116
SGM r = 5 0.00014 79 3.0 3.56 3.17 13-03 0.109
SGM r = 8 0.0009 96 4.2 7.29 3.12 13-03 0.097
SGM r = 10 0.001 102 4.9 9.85 3.08 E—03 0.101

 

 

 

 

 

 

 

 

 

Sinusoidal Flux (0' = (10°C)

 

 

 

 

 

GM whole 0.001 6 - 1.85 1.71 E-03 0.043
SGM r = 5 0.00012 80 3.1 3.48 1.70 13-03 0.044
SGM r = 8 0.0009 84 3.7 6.38 1.71 E-03 0.051
SGM r = 10 0.001 86 4| 8.34 1.71 E-03 0.050

 

 

 

 

 

 

 

 

 

 

Results from the analysis with exact data are given in Table 5-1 for estimating the
three heat flux functions. The first column identifies the method, whole domain gradient
method (GM) or a sequential gradient method (SGM). The second column lists the analy-
sis domain investigated, either whole domain or sequential with the value of r specified.
The Tikhonov regularization parameter is given in column three. Since exact data are

used, the regularization parameter is specified as 0.001 for the whole domain analysis and

118

varied to produce the same sum-of—squares error (column seven) for the sequential analy-
sis. The number of iterations, both total and average number per sequential interval, are
given in columns four and five. Computational time is listed in column six. Measures of
the error in temperature and heat flux are given in columns seven and eight, respectively.
The error in the temperature is computed as

1

J N 5
l 2
S}, = {mg 2 [T(dj,t,.)—Yj(ti)] } (5-4.2)

j=1i=l

Equation (5—4.2) represents a measure of the error in the sum-of-squares function J S in

equation (4-2.4). The residual (discrepancy) principle states that the error, S y in equation

(5-4.2), should be reduced to its expected value, Alifanov (1994). Error in the estimated
heat flux is

1

P N 5
A 1 A
CD = {mkx z [(1070 t,-)-61(rk, 1,112} (5-4.3)

= 11': 1

where q(rk, ti) is the true heat flux at location rk on F4 and time ti and 21 is the estimated
value with no measurement errors. Equation (5-4.3) is referred to as the deterministic bias,
or just bias, (Beck et al., 1985), since the estimated heat flux is computed without mea-
surement errors. The deterministic bias represents the bias that is introduced by the inverse
method itself; some bias is required to stabilize the ill-posed problem.

The estimated heat flux can more closely approximate the triangular shape for small
r-values using the SGM, particularly near the peak, see Figure 5—4. A smaller Tikhonov
parameter permits the SGM to more accurately reproduce the peak, with less bias or

smoothing of the estimated heat flux. It is interesting that the Tikhonov parameter is

119
reduced for the SGM at small values of r. The problem is more ill-posed for the SGM

compared to GM. Hence, it was anticipated that the Tikhonov parameter would increase,
but it decreases instead. The decrease in the Tikhonov parameter is required, even though
the problem is more ill—posed, to permit the solution the “flexibility” to reduce the sum-of-
squares to the magnitude specified by the residual principle. Regularization, tends to
“stiffen” the inverse problem. When solved in a sequential manner, estimates near the end
of the sequential interval are biased low (discussed below). To compensate for the biased
estimates near the end of the interval the Tikhonov parameter must be reduced to permit
“flexibility” in the first components and obtain the correct magnitude for the sum-of-
squares. Note that as r is increased, the Tikhonov parameter also increases towards the
magnitude specified for the whole domain solution.

Several observations are drawn from the analysis with exact data. The computational
time required for the sequential implementation was significantly more (by a factor or 2 or
greater) than the whole domain solution. Larger magnitudes of r required additional com-

putational time. Five was the minimum number of future times that could be used, which

represents a dimensionless time based on the sensor depth of t: = 0.3. Fewer than five

future steps, or more importantly a dimensionless time of 0.3, results in the estimated heat
flux being biased low. The additional computational time required and minimum number
of future time steps for the sequential implementation are explained by examining the
characteristics of the gradient method near the end of the time interval.

For illustrative purposes, consider that the steepest descent method is the iterative

search method used. In this case the search direction is (the negative of) the gradient,

which is related to the adjoint function. The heat flux at iteration n + 1 is

120

q" +1(ht) = q"(r, t) — 15571:, t). r 3 F4 (5-4.4)

For the steepest descent method this becomes

q“ 10‘, t) = qn(f, t)-(3"VJ(r, t),r3F4 (5-4.5)

Substituting for the gradient from the relationship developed for the adjoint function, gives

q" + 1(r, z) = q”(r, t) — p"{w(r, t) + aT[q"(r, t) — qpn-(r, 0]} , r 3 r4 (5-4.6)

Recalling that the adjoint function is prescribed to be zero at the final time, equation (4-
6.21d), equation (5-4.6) shows that updates to the heat flux approach zero at the final time.

Consequently, the estimated heat flux near the final time is biased. If no prior information

is used, q pn-(r, t) = 0, the estimated heat flux at the final time is approximately (the dif-

ference is the Tikhonov regularization term) the same as the initial guess (n = 0),

q" + 1(r, tm + r__ 1) z q0(r, tm + r _ 1) . For a zero initial guess the estimated value at the final

time is identically zero, q" + 1(r, t 0.

m + r— 1) =

It is not unreasonable to specify that the flux components near the final time approach
zero. Nor is it unreasonable to specify that the adjoint function be zero at the final time.
Very little information is available concerning the flux components near the final times.
Due to the lagging effect of the temperature response at the senor location to variations in
the surface heat flux, there is not time enough for information about the components near
the final times to reach the sensor. The magnitude of the lagging effect depends on the
dimensionless time (based on the sensor depth). For the same reason, the adjoint function,
which is shown to be closely related to the gradient, is insensitive to the components near
the end of the time interval. Appropriately, it is specified as zero there to simplify its com-

putation.

121

The inherent difficulty of estimating the heat flux near the end of the time domain is
more influential for the sequential implementation because a smaller number of time steps
are considered. As stated previously, for sequential intervals shorter than a dimensionless
time of 0.3 (five steps for these cases) results in biasing of the estimated heat flux. At the
threshold of 0.3, the response at the sensor and sensitivity is large enough such that the
estimated heat flux is not biased at the first time step of the sequential interval (the one
retained for a sequential interval). Components beyond the first time step are biased due to
the insensitivity of the temperature to these components.

An anticipated benefit of a sequential implementation is an improvement in the com-

putational efficiency. The expected improvement is based on the premise that after the first

sequential interval (t0 < t < tO + rAt) an accurate initial guess for the estimated heat flux

is available from the previous sequential interval. With an accurate initial guess it was
expected that few iterations, possibly one, would be required to converge after the first
sequential interval. Computational savings may not be realized for the one-dimensional
problem because the conjugate gradient method is quite efficient. Results indicate that
approximately the same number of iterations are required for all sequential intervals. A
result that indicates the conjugate gradient search method is insensitive to the initial esti-
mate. With a seemingly accurate initial estimate, convergence is expected within a couple
of iterations, especially for the one-dimensional problem and should improve after the first
sequential interval. This was not realized. the reason is discussed next.

The computational time required for the whole domain approach was 50-90% less
than the time required for a sequential implementation for r=5, see Table 5-1. As more

future time steps were considered. the computational time increases significantly, nearly

122

doubling for an increase from r=5 to r=8. Additional iterations, which in turn increases
the computational time, are required for the sequential implementation because of the
biasing in the estimated heat flux near the end of the sequential time interval. The biasing
of the estimates to zero (or a prescribed constant) results from the definition of the adjoint
function, which is prescribed to be zero at the final time. The sequential implementation
accentuates this effect, which was not expected to have such a significant effect. This was
not expected to be significant because only data at the first time step is checked for the
convergence. What was not realized was that in a sequential interval, the estimates at the
end of the time interval influence the estimates at the beginning of the time interval. Fur-
thermore, predicting the effect of the biasing is a difficult proposition. Its influence
depends on the variation of the (unknown) heat flux. Such an effect is difficult to predict or
anticipate.

Examining the estimates at each iteration for a sequential interval gives insight to this

effect. Estimated values for the triangular heat flux at each iteration for the sequential
interval beginning at I: = 0.12 are shown in Figure S-Sa, b, and c, for r-values of 5, 8 and

10, respectively. For r=5 in Figure 5-5a the initial estimate (n = 0) , set from values at the
previous sequential interval, is inaccurate. Biasing of the estimated heat flux resulted in a
40% error in this component during the previous sequential interval. At larger magnitudes
of r in Figure 5-5b and c, the effect at the beginning of the sequential interval is less pro-
nounced. However, the estimate at the first time step, which is the component retained for
a sequential interval, varies with iteration. For r=10 the converged estimate (at n=4) for
the first component varies only 2-3% compared to the initial estimate (n=0). The first com-

ponent varies 10-30% while iterating (n=1 and 2). Variations in the estimates are caused

123

 

 

 

 

 

0.80 ' I I I I I ' I ' I ' I ' I
G—ano
13—E1n=1
<>—<>n=2

3;, 0.60 r 1 “1:3
'QE’ — Exact
.2
a
0.)
.§.
3 0.40 — ~
é
a:
g 1
I
0.20 r .
0”8.00 I 0.10 L 0.20 ‘ 0.30 ‘040 I 0.50 l 0.60 I 0.701 0.80

+
t

Figure S-Sa Estimated one-dimensional triangular heat flux with no measurement
errors for sequential interval beginning at time 0.12 for r = 5

124

 

 

 

0.80
i3
.5
§
.§
3 0.40
§
Q:
§
I
0.20-
0“8.00 010 I 0.201 0.304 0.40 l 1
+
t

 

G—ano
[El—El n=1
<>—<>n=2
A—An=3
ak—alenz4
— Exact .

 

A;

0.50 0.70 L 0.80

Figure 5-5b Estimated one-dimensional triangular heat flux with no measurement

errors for sequential interval beginning at time 0.12 for r = 8

125

0.80'1e111r1vrr

 

0.60 .
E
g 0.40 -
§
.5 1
3
’5‘
E 0.20 ~
8
I

 

4 J 4

 

G—9n=0
El—Eln=1
<>—<>n=2
A—An=3 l
Hn=4
— Exact .

 

 

0.08

+
t

0010101020 0.3 0.4010.5010.601.7010.80

Figure 5-5c Estimated one-dimensional triangular heat flux with no measurement

errors for sequential interval beginning at time 0.12 for r = 10

126

by the biasing at the end of the time region influencing estimates throughout the entire
time region. In particular, estimates four to five time steps from the end (within a dimen-
sionless time of 0.3) have the largest change, and hence the most influence on the required
number of iterations. These components, not the one at the end of the interval, are the dif-
ficult ones. There is not enough sensitivity to estimate these component accurately, but
there is enough to cause them to vary. Similar results were seen when estimating a con-
stant heat flux.

Estimating the heat flux with exact data has indicated that a sequential implementa-
tion gives estimates within 5% (relative to maximum flux) of the whole domain method.
Computational time required for the sequential solution was significantly greater. Compu-
tational time is 50-90% greater than the whole domain, and increases depending on the
magnitude of r. Additional computational time is required because the sequential imple-
mentation is plagued by the inherent biasing of the components near the end of the time
interval. This issue is addressed later by introducing some functional forms (function

specification) over the sequential interval to control this effect.

5-4.2 Corrupted Data (Measurement Errors)

An important issue when investigating methods for solving the IHCP is the effect of
measurement errors on the solution. Methods that are typically insensitive to measurement
errors do so at the cost of introducing bias in the estimates. Conversely, methods that are
unbiased or have a small bias, are typically sensitive to measurement errors. These con-
flicting objectives, minimal bias and sensitivity to measurement errors, are further dis-

cussed later. The one-dimensional simulated test cases were evaluated again after adding

127

random measurement errors to the experimental data. Only errors in the measured temper-
ature are investigated. Errors have standard statistical assumptions and are described in
equation (5-4. 1).

Estimation of the triangular, step, and sinusoidal heat flux with measurement errors

added to the data are shown in Figure 5-6, 5-7, and 5-8, respectively. Results from the

 

 

 

 

 

0.8 v.vivr...,.
0——oGM
B—BSGMrzs
06~ HSGMr=8 ‘
' HSGM r=10
—Exact
.5 0.4— .
2
<1)
.§
3
E?
if 0.2~ -
S
I
0.0- -
_0.2 1 1 . 1 L 1 L 1 1 1 L
-1.0 -O.5 0.0 0.5 1.0 1.5 2.0
+
t

Figure 5-6 Estimated heat flux for one-dimensional triangular test case with
measurement errors

128

 

   
  

 

 

 

1.5 ' Tl ' I ' I I I ' I
e—eGM
9—9 SGM r=5
G—OSGM r=8
1.0 ~ A—ASGM r=10 ‘
:3 -— Exact
8
'E
.2
:5; 0.5 ~ -
é
C1:
6
m
0.0 - .
_05 . 1 L 1 . 1 L 1 L 1 L
-0.6 -0.1 0.4 0.9 1.4 1.9 2.5
+
t

Figure 5-7 Estimated heat flux for one-dimensional step test case with measurement
errors

129

 

  
  

Heat flux (dimensionless)

 

 

 

1.5 . . . , . 1 , r L 1
o—oGM
H SGM r=5
HSGM r=8
1.0 ~ 111—1186M r=10 ‘
— Exact
0.5 - -
0.0 — -
'0.5 1 1 1 1 L 1 1 1 1 1 1
-O.6 -O.1 0.4 0.9 1.4 1.9 2.5
+
t

Figure 5-8 Estimated heat flux for one-dimensional sinusoidal test case with
measurement errors

130

Table 5-2 Estimation results for simulated data corrupted with errors

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Iterations ,1
Analysis Tikhonov per (I32? SY 05'
Method domain (1.7 total seq int (sec) (°C) W/m2
Triangular Flux (0' = 0.0036°C )

GM whole 0.0028 6 - 1.84 3.63 13-03 0.034
SGM r = 5 0.0021 77 2.9 3.48 3.65 E—03 0.050
SGM r = 8 0.0019 91 4.0 6.95 3.64 E-03 0.040
SGM r = 10 0.0026 88 4.2 8.50 3.62 E-03 0.039

Step Flux (0' = 0.012°C)

GM whole 0.0033 7 - 2.06 1.23 E-02 0.139
SGM r = 5 0.0021 76 2.9 3.43 1.25 E-02 0.147
SGM r = 8 0.0026 87 3.8 6.66 1.24 E-02 0.1 19
SGM r = 10 0.0030 82 3.9 7.74 1.23 13-02 0.101

Sinusoidal Flux ( 0' = 0.006°C )

GM whole 0.0019 7 - 2.15 6.19 13-03 0.061
SGM r = 5 0.0016 76 2.9 3.42 6.19 13-03 0.088
SGM r = 8 0.0015 91 4.0 6.91 6.19 13-03 0.071
SGM r = 10 0.0019 92 4.4 8.88 6.19 13-03 0.070

 

 

 

 

 

 

 

 

 

analysis are given in Table 5-2. The tabulated results are similar to that presented for the
analysis with exact data. The type of analysis, whole domain (GM) or sequential (SGM),
Tikhonov parameter, number of iterations, computational time, and temperature and heat
flux errors are given in the table. In this case, which has measurement errors, the error in
the estimated heat flux represents an error due to bias in the algorithm and the influence of

measurement errors. This error is defined as the mean-squared error. The computation of

131

this error is identical to computation of the previous error in heat flux; however, the heat
flux is estimated with random errors present in the measurements
1

P N 2
A 1 A
GS, : {P(N _1) 2 2140700)" qe(rk2 ti)]2} (5-4.7)

k = 11': 1
where 218 is the estimated heat flux with measurement errors in the temperature.

The residual principle was used to select the magnitude of the Tikhonov regulariza-
tion parameter. As stated in the principle, the sum-of-squares function should be reduced
to the expected level. Since random errors are used, this is specified as the standard devia-
tion of the measurement errors.

The analysis with measurement errors, as expected, demonstrates similar computa-
tional aspects as the analysis without measurement errors, requiring significantly more
time for the sequential method. It also showed that the whole domain method was superior
to the sequential implementation in the accuracy of the estimated heat flux. Errors in the
estimated flux were up to 50% less, depending on the magnitude of r, and were largest for
smaller magnitudes of r. (Error in the estimated step heat flux was less for r—values of eight
and ten; however, this did not include the step down in the heat flux because this region
was within r steps of the final time.) These results show that, similar to the function speci-

fication method, the sequential gradient method requires a pr0per selection of r. However,
it further shows that selecting a magnitude of r and the Tikhonov parameter (r, (1T) , such
that the residual principle is met, does not insure optimum results. Improved results may
be obtained with another set of parameters. For smaller magnitudes of r the sequential

method is more sensitive to measurement errors, but the method reduces (IT to meet the

132

residual principle. Although the whole domain method provides a more accurate estimate,
the sequential method is competitive for a proper selection of r. If r is made large enough,
the estimates are independent of r and are the same for both the sequential and whole
domain solutions.

Quantifying the results of an estimation are discussed in Beck et al., (1985) and Scott

and Beck (1989). The errors suggested for quantifying the results of an estimation are the

deterministic error D and the mean-squared errorSe. These errors provide insight into the

characteristics of an estimation algorithm. The errors are computed in equation (5-4.3) and
(5-4.7) and are tabulated in Tables 5-1 and 5-2. For the deterministic error the heat flux is
estimated with no measurement errors, while for the mean-squared error the heat flux is
estimated with random errors in the measurements.

The deterministic errors, which show the bias of an algorithm, are smaller for the
sequential approach for the r=5. The deterministic bias, in the final column of Table 5-1, is

smaller because OtT is smaller for all three cases for r=5. Since there is less regularization

the solution can more quickly respond to changes in the heat flux. However, as r increases

(and a7 increases) the deterministic bias for the SGM becomes larger, and equals or

exceeds the whole domain value. The mean-squared error, which can be compared to the
deterministic error to determine the sensitivity to measurement errors, is exclusively better
for a whole domain solution. The mean-squared error is shown in the final column of
Table 5-2. A sequential method permits the Tikhonov parameter to be smaller at low r—val-
ues, however, this is at the cost of an increased sensitivity to measurement errors, as shown
by a larger mean-squared error at these r—values. As r increases the SGM becomes less

sensitive to measurement errors, of course in the limit the SGM and GM are the same.

133

Although the SGM provides comparable results to the GM, the computational
requirements are significantly greater. A 100% (or more) increase in the computational
time is typical for a sequential implementation compared to a whole domain solution. The
increase in computational time is a direct consequence of the sequential implementation
(of length r—steps) requiring approximately half as many iterations for each sequential
interval as the whole domain requires for a complete solution. Since these results are for
the one—dimensional problem, when the methods are extended to two-dimensions the
requirements for the sequential implementation might be greater in comparison to the
whole domain approach. As demonstrated previously, the additional computational
requirements are caused by the inability to estimate components near the end of the
sequential interval.

Several approach were investigated to improve the computational requirements of the
sequential implementation. It was clear from observing the results during a sequential
interval (see Figure 5-5) that the variation in components near the end of the interval result
in an increased number of iterations. If these components over the end of the sequential
interval could be controlled, or anticipated in their response, improvement in computa-
tional requirements is possible. Note that the components in the range 4-5 steps from the
end of the interval (a dimensionless time of 0.25 to 0.3) are troublesome. On either side of
these components the response of the heat flux is understood.

Approaches to improve the computational efficiency of the sequential method looked
at three aspects of the sequential implementation. First is the setting of the initial condition
for a sequential interval. Next, was using regularization to control the variation of certain

components in a sequential interval. Lastly, was the use of prescribed functions over the

134

sequential interval, similar to the function specification method. As it turned out, only the
final approach, using function specification over the sequential interval significantly
improved the computational efficiency. Much was learned concerning the method, how-
ever, and is described below.

Setting the initial condition for the subsequent sequential interval, from the previous
interval, has all the required information except for I-components, where I is the number
of components retained on time for a sequential interval. Initially, all values from the past
interval were carried to the next interval then the final I-values were specified to be zero.
But it was noticed that the values near the end of the interval did not change with sequen-
tial interval. Another procedure was used. Estimates from the previous interval were used

from the beginning and end of the sequential interval. The procedure is shown graphically

in Figure 5-9. For a specified number of time steps, m end (typically m en d = 4) from the

end of the interval, the initial values were set from the previous interval

0
q (r’t) : qNC0"V(f,I—IAI) ’tm+r-1—mmd+lStStm-+r—l (5-4'8)

where qum’ is the heat flux estimated at the previous sequential interval and q0 is the
initial estimate for the next interval. Equation (5-4.8) specifies that the initial values at the
end of the sequential interval should be the same as they were at the end of the previous
sequential interval. The components are shifted I-values in time because the time interval
advances l—values for the next sequential interval. Initial values at the beginning of the

time interval are specified directly from the previous interval

0
q (r, I) = qNCOflV(’-,t) ’tm—1+IStStm+r-l—m (5-4.9)

rnd

135

1!

tm+r—l+l

 

 

 

N
/ I.
N / I
/ + 1.
// / E
_1 /
j / / g
1. / /
l’ +/ /
/ E
/ s
‘7 / 1
:,/ .l
+' ‘ I
E L f
“ '- /’ 1 .5:
1 /
I '/ T a
. /'
.3 // .1 o
E , ..1
E // - 33
l 1
_ / ---- 45
| / ' t:
1- " .1 ' o
+_ + . j‘ S1
E 1 ,2
N 1 E” o
3 '1 1 U
\1 . .... I
E .
« l _______ #L 1'. “N
‘~ +
+
5 NE
N
--~
I -1
:_ -+
H)- + —————— h;- “E
E
H
‘
-1~ E ,...
‘

 

 

U"
(1 'SIMIOJNb

Figure 5-9 Procedure for specifying the initial condition using converged estimates
from the previous sequential interval

136

This leaves I-values (where I is the number of retained components on time) to be speci-
fied, which are computed using linear interpolation to join the two previously defined
regions.

Unfortunately, accurately specifying the initial condition from the previous interval
was not causing the increased computations. Computational requirements were unchanged
using the procedure outlined for setting the initial estimated. Rather, it was understanding

that the estimated heat flux in the previous sequential interval was not as accurate as antic-

ipated, particularly components about 4-5 steps (or t: = 0.3) from the end of a sequential

interval. These components vary significantly from one sequential interval to the next,
which influences other components in the interval.

The results in Figure S-Sa, b, and c show that the initial condition was not as accurate
as believed. Furthermore, the components that cause the most difficulty in the sequential
estimation are the last 4-5 components (covering the final 0.3 in dimensionless time of the
interval). Consequently, regularization was used to attempt to control the change in these
components. The Tikhonov regularization parameter was specified to be a function of
time, with magnitudes larger near the end of the sequential interval. The larger regulariza-
tion parameter drives the estimated heat flux to zero. This approach did eliminate changes
in the component near the end of the sequential. But, the changes that would have
occurred in these component migrated to other components. In retrospect this should have
been anticipated, since the solution for the heat flux should conserve the total energy (for a
specified r-value), by forcing some components to be zero (or a constant), other compo-
nents will accommodate to conserve the total energy. Computational aspects were not

improved.

137

A positive aspect of the gradient methods is their function estimation basis. The diffi-
culties encountered in the sequential implementation, however, are due to the attempt to
estimate a function over the sequential interval, even though, only one component of the
function is typically retained. A final approach to reduce the computational requirements
is to constrain the heat flux over the sequential interval so that not as much computational
effort is needed to estimated all components on a sequential interval. Fewer parameters are
estimated instead of a function. Specifying several functional forms over the sequential
interval were tried: cubic, parabolic, linear, and constant. Additional constraints were
imposed for the higher order functions, such as zero slope or zero second derivative at the
end of the sequential interval, as required. Prescribing a single function over the sequential
interval did not improve the computational aspects when higher order functions were used.
However, specifying the heat flux be constant over the sequential interval did improve
computational efficiency significantly.

The triangular heat flux case was re-evaluated using function specification over the
sequential interval. Herein, this method is referred to as the sequential function specifica-
tion gradient method (SFSGM). Maintaining the heat flux constant over the sequential
interval stabilizes the problem and the Tikhonov regularization parameter is set to zero.

The results for the exact data are shown in Table 5-3 for the triangular heat flux using
the function specification assumption; graphical results are not shown. The accuracy of the
estimated heat flux is comparable to the whole domain. A direct comparison is not made
because the Residual principle is not met by all r-values. The magnitude of the sum-of-
squares error suggest an r between 4 and 5. In this range the error in the estimated heat

flux is comparable, and more importantly computational requirements are reduced. For

138

Table 5-3 Estimation results for prescribing a constant heat flux over the sequential

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

interval
Analysis Tikhonov 1“”er C113: SY GD 0" 63'
Method domain aT total seq int (sec) ( o C) W/m 2
Triangular Flux (G = 00°C)

GM whole 0.001 7 - 2.14 8.59 E-04 0.0094
SFSGM r = 3 0.0 51 1.8 1.08 9.65 E-05 0.0054
SFSGM = 4 0.0 49 1.8 1.47 6.07 E-04 0.0052
SFSGM r = 5 0.0 48 1.9 1.90 1.89 E-03 0.014
SFSGM r = 6 0.0 48 1.9 2.24 9.13 E-04 0.025

Triangular Flux (0' = 0.0036°C)

GM whole 0.0028 6 - 1.84 3.63 E-03 0.034
SFSGM = 4 0.0 54 2.0 1.62 2.52 E-03 0.096
SFSGM r = 5 0.0 50 1.9 1.98 3.27 E-03 0.039
SFSGM r = 6 0.0 50 2.0 2.38 4.99 E-03 0.040
SFSGM r = 8 0.0 46 2.0 2.93 1.07 E-02 0.048

 

 

 

 

 

 

 

 

 

 

r=5, using the function specification requirement (SFSGM) the computational time is
reduced by nearly 50% compared to the standard SGM; computational time is 10% less
than a whole domain solution.

The heat flux estimated with measurement errors added to the data is shown in Figure
5-10 when assuming the heat flux is constant over the sequential interval. Estimation
results are also given in Table 5-3. With the heat flux held constant over the sequential
interval the sequential implementation improves significantly in computational efficiency.

Since only one components is estimated at each sequential interval, convergence is

139

obtained with few iterations. When the heat flux is held constant over the sequential inter-

val, denoted SFSGM, a sequential solution is competitive with the whole domain solution.

 

 

 

 

 

0.8 ' I ' l u T ' I 1 I r
o—oGM ‘
H SFSGM r=4
0 6 _ " <+—<>SFSGM r=5
”a? ' A——ASFSGM r=6
é’ HSFSGM r=8
.g — Exact
t:
E
g 0 4 b ‘
><
=2
5
m 0.2 - -
I};
.\
0.0 ” L‘3'16'.5; : 3 -
‘.
_02 1 1 1 1 1 1 1 1 L 1 1
-1.0 -0.5 0.0 0.5 1.0 1.5 2.
+
t

Figure 5-10 Estimated heat flux for one-dimensional triangular case with the heat flux
held constant over future sequential interval

140

The computational time is further reduced for a sequential solution by using prior
information. The prior information enters the solution in the Tikhonov regularization term
of equation (4-2.4). Specifying the prior information as the initial estimate, which is set
from values at the previous sequential interval, is shown to significantly improve the accu-
racy of the estimated heat flux as well as reduce the computational time required for a
sequential solution. The use of the prior information is discussed in the next chapter.

Results using prior information are not shown for the one-dimensional case.

5-4.3 Results - Experimental Measurements

Although simulated measurements provide insight to the behavior of the method, the
true test is with experimental data. Gradient and iterative regularization methods are
applied extensively in the literature to solve the IHCP. However, few of these investiga-
tions have applied the method to experimental data. Measured temperature and heat flux
histories from the experiment to estimate the thermal conductivity are used. For this pur-
pose, the problem is reversed and the estimated thermal properties and measured tempera-
tures are used to recover the transient surface heat flux. Since electric heaters were used to
generate the heating, the heat flux is also quantified. Consequently, a comparison between
measured and estimated heat flux can be made.

The experimental apparatus is shown in Figure 2-1 and discussed in section 2-2.0.
Electric heaters symmetrically heat the composite material while thermocouples measure
the transient temperature. Since the experiment is symmetric, the measured power to the
heaters quantifies the energy and heat flux to the specimens. In the thermal model of the

experiment, shown in Figure 2-2, the thermal effects of the mica heater, and more

141

importantly the contact resistance, are included. An effective thermal conductivity, which
represents the heater, contact resistance, and the cement used to install the thermocouples,
is used to group these effects. Heat losses to the insulation are also included. Temperature
measurements are made very close to the surface where the heat flux is estimated (a depth
of 0.44 mm). However, the relatively low effective thermal conductivity of the mica
heater/contact resistance results in the problem being more difficult than it appears. The
Fourier number based on the depth of the closest sensor and thermal properties of the mica
heater/contact resistance is approximately 0.06. The relative high thermal conductivity of
the carbon-carbon specimen increases the difficulty of estimating the surface heat flux
also. This is demonstrated by the limiting case, if the conductivity of the specimen were
infinite, the sensitivity to the surface heat flux is zero and it is impossible to estimate it.
This case was also investigated by Beck et al. (1996), which among other things
investigated the use of the residual principle for selecting the number of future time steps
for the function specification method. Since seven reliable thermocouples were used to
measure the temperature at nominally the same location, the error in the measurements
could be quantified. Beck et al. suggested using the estimated standard deviation of seven

measurements from the average temperature at each time. The standard deviation was
approximately 004°C at the beginning and end of the experiment when the heater was

not active and 012°C when the heat flux was active. A time average of this standard devi-
ation was 0.083 C. This value was used as the level of noise for selecting the regularization
parameters with the residual principle.

Beck et al. (1996) also compare several inverse methods. Function specification,

Tikhonov regularization, iterative regularization, and specified functions over large time

142

regions using Green’s functions were compared. Results indicate that estimates from the
various methods are quite comparable. Although computational requirements were not the
focus of the study in Beck et al. (1996), it was noted that certain methods require consider-
ably more computational time.

Estimations of the surface heat flux are shown in Figure 5-11 for the experimental
case. Several variations of the gradient method are compared to estimate the surface heat
flux. Whole domain gradient method and sequential gradient method including specifying
a constant flux over sequential internal are shown. A remarkable outcome of the different
techniques is that about the same results is obtained. Results of the analysis are also listed
in Table 54. Only for r=4 is there a significant difference between the whole domain and
sequential solutions. Otherwise, the sequential solutions are as good as, or better than, a
whole domain solution. The sequential solutions, denoted SGM in column one of Table 5-
4 require about three iterations per sequential interval. Holding the heat flux constant over
the sequential interval, denoted SFSGM in column one, requires about two iterations per
sequential interval. This represents a reduction of approximately 30% in computations
with SFSGM. A larger decrease in the computational time is achieved because the number
of iterations required for SFSGM is actually less. Although the second iteration is started,
frequently the gradient computed during the second iteration is very small (less than 1E-
15) and based on the magnitude of the gradient iteration is stopped. So a full iteration is
not completed (the sensitivity problem is not solved) and greater savings are realized.

These additional savings increase with r.

143

 

 

 

 

 

 

 

 

 

 

 

 

 

9000.0 1 1 1 1 1
o—oGM
B—EISGM r=8
— Measured
A5000.0 - -
”E 1
E
6
{13 3000.0 - , -
8
I J
1000.0 - .
-1000.0 1 ‘ 1 1 ' '
0.0 10.0 20.0 30.0 40.0

time (seconds)

Figure 5-11 Estimate heat flux for one-dimensional experimental case 1010#30.l

144

Table 5-4 Estimation results for analysis of experimental case

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Iterations
Analysis Tikhonov per (£1912? SY 64
Method domain a7" total seq int (sec) (0 C) W/m2
Experiment Case 1010#30.1
GM whole 1.0 E-08 8 - 24.1 0.080 766
SGM r=4 4.0 E- 10 603 3.1 30.0 0.079 1400
SGM r=6 2.0 E-09 564 2.9 47.4 0.077 754
SGM r=8 4.5 E-09 560 2.9 64.7 0.079 716
SGM r=10 7.0 E-09 565 2.9 80.2 0.078 726
SFSGM 124 0.0 382 1.9 16.5 0.039 2250
SFSGM r=6 0.0 373 1.9 26.4 0.044 762
SFSGM r=8 0.0 347 1.8 33.1 0.057 699
SFSGM r=10 0.0 320 1.7 39.4 0.083 751
SFSGM r=6 3.5 E-lO 370 1.9 25.4 0.086 764
SFSGM m8 4.5 E-10 342 1.8 32.7 0.077 714
Experimental Case 1023 #4031
GM whole 4.0 E-09 8 - 25.3 0.080 1230
SGM i=6 9.0 E- 10 578 3.0 47.0 0.084 1182
SGM :8 2.0 E-09 578 3.0 82.1 0.084 1141
SGM :10 3.0 E-09 576 3.0 82.5 0.079 1167
SFSGM r=6 0.0 359 1.9 26.0 0.036 1234
SFSGM r=8 0.0 342 1.8 34.2 0.052 1162
SFSGM r=10 0.0 306 1.8 39.8 0.081 1261
SFSGM r=6 1.5 E- 10 359 1.9 26.0 0.080 1215
SFSGM r=8 2.0 E- 10 336 1.8 32.6 0.077 1180
SFSGM r=10 1.2 E- 10 307 1.6 39.9 0.084 1267

 

 

 

 

 

 

 

 

 

Chapter 6

APPLICATION OF THE SEQUENTIAL GRADIENT METHOD:
TWO-DIMENSIONAL IHCP

6-1.0 Introduction

The solution of the two-dimensional IHCP is investigated in this chapter. Use of the
sequential gradient method (SGM), possibly coupled with prescribing a functional form
over the sequential interval (SFSGM) are studied. A comparison of sequential gradient
methods with the standard whole domain gradient method (GM) is made. Some compari-
son with the standard function specification method, which does not use a gradient search
method nor an adjoint equation approach, is presented. The two (or three) dimensional
problem is the intended use of the proposed sequential methods. As demonstrated with the
one-dimensional cases, assuming the heat flux is constant over the sequential interval
improves the computational time of the sequential implementation, and makes it competi-
tive with the whole domain approach. Additional steps may be required to improve the
computational aspects for the two-dimensional problem also.

An approach similar to the one-dimensional analysis is pursued for the two-dimen-
sional case. Several test cases for which the surface heat flux is known are studied to quan-
tify the accuracy of the method for a two-dimensional analysis. The two-dimensional

cases will demonstrate the ability to estimate spatially dependent fluxes. In addition to

145

146

simulated cases, an experimental case with a measured surface heat flux, which is two-
dimensional and transient, is analyzed.

Several variations of the sequential gradient method are suggested in this chapter.
Three methods are investigated. The standard whole domain gradient method (GM) solves
for the heat flux with data from all measurement times. The sequential gradient method
(SGM) solves for the heat flux using data from a sequential interval. The sequential func-
tion specification gradient method (SFSGM) solves for the heat flux over a sequential
interval while prescribing a functional form for the heat flux on time, but allows for spatial
variation. In addition, the use of prior information is suggested during the sequential solu-
tion to improve the sequential method. This is the first known application of prior informa-
tion to improve a sequential solution.

The remainder of this chapter is summarized. The calculation of the simulated tem-
perature data is discussed in the next section. Section 6-3.1 presents results for the simu-
lated test cases with exact data. Test cases with measurement errors added to the simulated
data are discussed in Section 6-3.2. An experimental case, with a measured surface heat

flux is analyzed in Section 6-4.0.

6-2.0 Simulated Temperature Data

A rectangular domain of dimensions [0, a] x [0, b] is studied for the two-dimen-
sional simulated cases. A schematic of the geometry is shown in Figure 6-1. One bound-
ing surface of the two-dimensional body is assumed to have an unknown heat flux while
all other surfaces are insulated. Furthermore, constant anisotropic thermal properties and a

uniform initial temperature are assumed.

147

q(x1t)

 

 

 

 

/‘ \ ’
4 O O O O O O O O Jk T

o / \ o

11 6 § 11 b
5 x h
/ \\\\\\ \\\\\\ \\\\\\ \\\\\\\\\\\\ \\\\\\ \\\\\ .\ I

q = 0
y f a $»

 

 

 

 

Figure 6-1 Two-dimensional geometry for simulated test cases

The problem can be mathematically written as

 

2 2
[bzkaa T: a 7;" 6T?“
+ — _

 

__ _ (6-2.10)
aZky ax+2 ay+2 81+
aT?‘
_ a i = (13,114", 1*) = ? (6-2.lla)
.v ),.:0
8T:
= 0 (6-2.llb)

 

+
3y

148

arf
; = 0 (6—2.l 1c)
x
x+ = 0
Br:
I
x+ = 1
T:(x+, y+, 0) = O (6-2.lle)

where the following dimensionless variables were introduced

 

q: = {-15}; (62.12)
T-—T
+ l 0
= —— (6-2.l3)
' (qNb/k),)
k,/ c)!
+ = ( " z (6-2.l4)
b
x+ = (1: (6-2.15)
y+ = 11) (6-2.16)

The subscript 1' indicates the prescribed heat flux. The estimation of two different pre-

scribed heat fluxes is investigated for the two—dimensional case. A step change in heat flux

on time (starting at t+ = 0 and ending t+ = t2) and position (from x+ = 0 to x)r = a1)

is the first flux (1' = l) examined

qi(x+,t*)= l (O<t+<t;),(OSx+Sal) (64.17)

0 otherwise
The second flux ((i = 2) examined also has a step change in time, but the spatial profile is

triangular over the entire surface at (y = 0)

149

213 (0 < 1" < 1;),(0 s x“ s

l

2—211+ (0<t+<t;),G<x+S 1)

Nil—

q:(x+, t+) = < (6-2.18)

 

1 0 otherwise

Due to the complexity of obtaining analytical solutions for the two-dimensional prob-
lems, an accurate finite element code (TOPAZ2D) is used to solve for the two-dimensional
temperatures. Data for the simulated cases are shown in Figure 6-2 and 6-3 for the step
heat flux and triangular heat flux, respectively. Temperature data at the sensor locations

and the heat flux at the surface are shown for the two cases. Both surface flux histories

\

l
begin at time 0 and end at time t;, which is prior to the end of the time domain, t}. The

temperature data are generated for eleven simulated sensor locations that are equally

spaced along the x-plane at a depth of y+ = 0.1 . Lagging effects are similar to the one-
dimensional problem. Notice that the spatial temperature profile is smoothed compared to
the profile of the surface heat flux. The diffusive nature of the problem at the sensor depth

makes it difficult to recover sharp spatial features.

6-3.0 2D Results - Simulated Measurements

6-3.l Exact Data (No Measurement Errors)

A case with “exact” data is analyzed first. The dimensionless groups of interest for

this geometry are a/b = 2 and kx/ky = 1 . Measurements are available at eleven equally

spaced locations (x+ = 0.0, 0.1. . . ..l.0) beneath the surface of unknown heat flux at a

depth of y+ = 0.1. The dimensionless time step (Fourier number) based on the sensor

150

fig

LLLLLL 2%

//
__LLLL///////////////////////W//////////MM2

/
2222/2222

3::

//

//
222/2

//

//

2222

 

\
2
0

AmmoEommcoEEV BBEQQEDH

////////o

 

 

\ \ \ \ w >
1 8. 5 A 2 05
o 0 0 0 1

Ammo—coficoEé 5E “mom oomtzm

 

Figure 6-2 Simulated temperature and surface heat flux for a two-dimensional

triangular heat flux

151

0.2

:2
:g

:

,
\LLLL_____:%
:52???
£22

55:2

,
:3:

:,

     

  

—5

  

/
/
7;

 

0.3 \

0.2 \

Ammo—coicoEEv 11:58an .H.

 

 

 

 

Ammo—coicoEEV xsc Hao: oomtsm

 

Figure 6-3 Simulated temperature and surface heat flux for a two-dimensional step heat

flux

152

depth is At: = 0.06 , which represents a difficult case; two-hundred time steps are consid-

ered in the analysis. “Exact” simulated temperature data are generated from the finite ele-
ment solution using TOPAZZD (Shapiro, 1986)

The geometry shown in Figure 6-1 is discretized for the numerical solution. A finite
control volume mesh with eleven nodes in the x-direction and twenty-one nodes in the y-
direction is used. A more coarse node spacing is specified for the x-direction. Because
eleven sensors are located along this surface, eleven components of heat flux are to be esti-
mated. It is possible to define more nodes along this surface and estimate a greater number
of spatial components (discussed later). Altemately, more nodes can be defined along the
surface and fewer spatial components, compared to the number of nodes, can be estimated
by introducing some spatial smoothing using basis functions; see Section 4—8.0. In this
study, it is not desired to impose spatial smoothing on the solution. Instead the number of
nodes defined for the numerical solution are set based on the number of available sensors.
Hence, eleven nodes are defined on the surface. The number of nodes is increased and
more spatial components are estimated in a later analysis. Care is exercised to not have too
few of nodes and incur numerical error from a coarse discretization.

The estimated surface heat flux for a whole domain solution (GM) and the sequential
gradient method (SGM) with r=6 are shown in Figure 6-4a and b, respectively, for the tri-
angular test case. The GM has the whole domain characteristic of over—shooting near
regions where the heat flux has an abrupt change, at the beginning and ending times. The
SGM reduces the effect of these characteristics. Both methods accurately estimate the sur-

face heat flux.

Heat flux (dimensionless)

153

a) GM
Whole Domain

0.8 4

0.6 \1

0.4 x

0.2 \

 

 

 

 

 

 

15

 

 

 

0.8

 

 

Figure 6-4a Estimated surface two-dimensional heat flux for triangular test case using
data without measurement errors. Whole domain solution

Heat flux (dimensionless)

154

b) SGM
r = 6

 

 

 

Figure 6-4b Estimated surface two-dimensional heat flux for triangular test case using
data without measurement errors. Sequential solution with r=6

155

Table 6-1 Estimation results for the two-dimensional IHCP with exact simulated data

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Analysis
domain Iterations
Tikhonov CPmP S Y 6 D
per t1me
Method 1' I “T total seq int (sec) (°C) W/m2
Triangular Flux (11 Components)

GM all all 1.00 E04 1 l - 35.7 2.53 E-04 7.46 E-02
SGM 6 1 1.02 E-04 758 4.0 71.8 2.50 E-04 8.05 E-02
SGM 10 1 3.20 E04 735 4.0 121.2 2.52 E-04 7.79 E-02

SFSGM 6 l 1.05 E-05 244 1.3 23.2 2.51 E-04 7.89 E-02
SFSGM 10 l 0.0 238 1.3 35.6 4.29 E-04 7.32 E-02
Step Flux (11 Components)

GM w all 1.0 E-04 12 - 38.8 4.42 E-04 9.30 E-02
SGM 1 1.25 E-04 762 4.0 71.7 4.46 E-04 9.92 E-02
SGM 10 1 4.05 E-04 737 4.0 119.9 4.43 E-04 9.54 E-02

SFSGM 6 1 1.67 E-05 264 1.4 24.9 4.42 E-04 9.68 E-02
SFSGM 10 l 0.0 260 1.4 42.1 5.41 E-04 8.84 E-02
Triangular Flux (21 Components)

GM w all 1.00 E-04 11 - 66.3 2.53 E-04 7.46 E-02
SGM 6 1 1.01 E-04 758 3.9 130.7 2.50 E-04 8.12 E-02
SGM 10 1 3.10 E04 727 3.9 222.4 2.50 E04 7.72 E02

SFSGM 6 l 1.05 E-05 270 1.4 47.0 2.52 E-04 9.65 E-02
SFSGM 10 1 3.00 E-OS 278 1.5 83.6 4.91 E-04 7.89 E-02

 

 

Pertinent information about the estimations with exact data are given in Table 6-1.
The table lists the solution method and analysis domain, with the number of future time
steps and number of components retained on time in the first three columns. Column four

gives the magnitude of the Tikhonov parameter. The number of iterations, both total and

156

average per sequential interval are given in columns five and six. Computational time is
listed in column seven. The errors in the sum-of-squares temperature and estimated heat
flux (bias error) are shown in columns eight and nine. Values of r = 6 and r = 10 were
studied for the number of future time steps. The Tikhonov regularization parameter was
specified as 0.0001 for the whole domain solution (GM). Subsequent sequential solutions
(SGM and SFSGM) varied the Tikhonov parameter to match the magnitude of the sum-of-

squares, S y, to the level of the whole domain solution.

Estimation results with exact data for the triangular and step heat flux, given in Table
6-1, have trends similar to the one-dimensional case. A comparable heat flux to the whole
domain (GM) solution is estimated with the sequential implementation (SGM), the bias
error (last column of Table 6—1) is 6-8% larger for r = 6, but it requires approximately
twice as much computational time. Increasing the r-value to ten reduces the bias error to a

magnitude 2-4% larger than the GM, but requires a proportional increase in the computa-

tional time, i.e. computational time increases by 167%(10/6) , compared to the computa-

tional time for r = 6.

In the one—dimensional solution the computational aspects were improved by main-
taining the heat flux constant (on time) over the future sequential interval, which is the
standard function specification assumption. The same procedure was tried for the two-
dimensional problem. Prescribing that the heat flux remain constant over the sequential

interval (method SFSGM in Table 6-1), requires one-third the computational time of the
SGM for r = 6, and has a bias error that is only 4-6% larger than the whole domain esti-

mation. Improved accuracy was obtained when r was increased to ten, bias error was 2-5%

157

less than whole domain solution. Notice that the Tikhonov parameter was set to zero for
the SFSGM with r = 10. Even with this parameter set to zero the sum-of-squares was not
reduced to the required level. This is due to the effects of the Tikhonov parameter and the

function specification approximation, which are discussed next. The SFSGM has superior
computational aspects. It requires 35% less computational time than the GM for both the
triangular and step test case for r = 6; computational requirements are comparable for
r = 10.

One observes in Table 6-1 that the estimation with exact data requires a larger magni-
tude for the Tikhonov parameter when solving sequentially compared to whole domain.
Since a larger magnitude of the Tikhonov parameter is used, the sequential solution is
more biased than the whole domain solution, particularly near regions of variation in the
heat flux. An increase in the Tikhonov parameter is expected because the shorter duration
sequential problem is more ill-posed than the whole domain problem. For the one-dimen-
sional solution, in general, a smaller Tikhonov parameter was required for the sequential
solution compared to the whole domain solution, contrary to what was expected. It is
believed that the one-dimensional problem required a smaller Tikhonov parameter to per-
mit the solution algorithm the flexibility to reduce the sum-of—squares magnitude to the
desired level. The reduction in the Tikhonov parameter is greater at small r-values. As r is
increased the Tikhonov parameter for SGM approaches the magnitude of GM, in the one-
dimensional case.

Applying the function specification approximation in the inverse solution provides
significant regularization, which increases with r. As noted in Chapter 5, regularization

imposes “stiffness” to the inverse solution. When maintaining the heat flux constant over

158

the sequential interval, considerable “stiffness” is imposed on the solution. Consequently,

the Tikhonov parameter must be reduced in magnitude to permit obtaining the desired

sum-of-squares magnitude, S y. At some level of r, the bias error becomes large enough

such that even when the Tikhonov parameter is specified as zero, the sum-of-squares func-
tion can not be reduced any further. This is a result of the bias introduced when assuming
the heat flux is constant over the sequential interval.

A noted benefit of the gradient method is that it is not required to specify the func-
tional form of the unknown heat flux. Since a numerical solution is used, however, some
approximation of the heat flux is required. The simplest numerical approximation of the
heat flux is used, which assumes the heat flux is constant over a control surface surround-
ing the node. Consequently, the functional heat flux is numerically represented by P spa-
tial components, where P is the number of nodes along the surface of the unknown heat
flux. Specifying more nodes along this surface will improve the numerical accuracy of the
direct problem, but requires estimating more spatial components, which may increase the
difficulty of the inverse problem depending on sensor locations. A general rule is that the
number of spatial components estimated should not exceed the number of sensors. Fewer
spatial components of the heat flux, compared to the number nodes, can be estimated by
prescribing a spatial functional form using basis functions, see Section 4-8.0. Prescribing
a spatial functional form introduces regularization and helps stabilize the solution.

In the previous simulated cases, given in the first two blocks of Table 6-1, eleven
nodes were defined along the x-direction and eleven spatial components were estimated.
Now the triangular test case is analyzed with twenty-one nodes components on the surface

and twenty-one spatial components are estimated, with eleven equally spaced sensors

159

along the surface. Results are given in the final block of Table 6-1. The estimated heat flux
is nominally the same for the whole domain (GM); of course computational time is
approximately doubled with twice as many nodes. The sequential method does not per-
form as well. Bias error is increased by 1% for the SGM, but increases 22% for the
SFSGM with r=6. When estimating more spatial components than available sensors with
a sequential method the spatial components between the sensor locations are biased low
for r=6. This produces spatial oscillations in the estimated heat flux, a “zig-zag” pattern
about the true flux value exists in the estimated heat flux. The oscillations are reduced by
increasing the magnitude of r.

It appears that estimating more spatial components than available sensors is possible
using the gradient method. At least for the case when sensors are evenly distributed along
the surface it is possible. The method seems to have some spatial regularization inherent.
Implementing a sequential solution degrades this characteristic for small values of r.
These outcomes are influenced by the dimensionless time step based on the sensor depth,
the spacing of the sensors along the surface relative to the sensor depth, and the amount of

noise in the measurements, which was exact for this case.

6-3.2 Corrupted Data (Measurement Errors)

The two-dimensional test cases are analyzed with measurement errors added to the
simulated temperatures. Standard statistical assumptions are used to describe the measure-
ment errors. The errors are modeled as additive, with zero mean, constant variance, and a

normal distribution. A magnitude of approximately 1% of the maximum temperature rise

160

was selected as the standard deviation of the measurement errors. This error corresponds
to 0.0018°C and 0.0025°C for the triangular and step heat flux, respectively.

Estimation of the triangular heat flux with measurement errors is shown in Figure 6-
5a for a whole domain solution and Figure 6-5b shows a sequential solution with r=6. The
two solutions are quite different, whereas the whole domain solution is relatively smooth,
the sequential solution has a significantly larger variability. It is not clear if the variability
of the SGM estimated heat flux is larger on time or on space. The sequential results are

certainly unacceptable in comparison to the whole domain solution for r=6.

a) GM
Whole Domain

1 \
0.8\

0.6 w

 

Heat flux (dimensionless)

 

 

 

 

 

 

 

Figure 6-5a Estimated surface heat flux for triangular test case using data corrupted
with measurement errors (0' = 0.0018°C ). Whole domain solution.

161

 

 

sing data corrupted

with measurement errors (0' = 0.0018°C ). Sequential solution with r=6.

Figure 6-5b Estimated surface heat flux for triangular test case u

162

Table 6-2 Estimation results for two-dimensional IHCP with simulated data corrupted
with random errors

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Analysis
domain Iterations A
Tikhonov Comp S y as,
per t1me
Method r I (1.7- total seq int (sec) (0 C) W/m2
Triangular Flux (0' = 0.0018°C) - 11 Spatial Components

GM all all 1.70 E-03 9 - 31.8 1.82 E—03 8.57 E-02

SGM 6 1 2.20 E-04 764 4.0 77.8 1.81 E-03 1.36 E-Ol

SGM 8 l 4.50 E-04 758 4.0 96.0 1.80 E-03 1.02 E-OI

SGM 10 1 7.00 E-04 746 4.0 1 19.0 1.80 E-03 9.19 E-02

SGM 15 1 1.25 E-03 718 4.0 173.7 1.80 E-03 8.65 E-02

SGM 20 1 1.50 E-03 759 4.3 246.9 1.80 E-03 7.53 E-02
SFSGM 6 1 5.00 E-05 767 4.0 77.1 1.81 E-03 1.57 E-Ol
SFSGM 8 1 8.30 E-05 769 4.1 96.3 1.80 E-03 1.12 E-01
SFSGM 10 1 1.1 E-04 761 4.1 120.2 1.80 E-03 9.71 E-02
SFSGM 15 l 9.50 E-05 693 3.8 167.2 1.87 E-03 8.97 E-02

Step Flux (0' = 0.0025 ° C) - 11 Spatial Components

GM all all 6.00 13-04 13 - 45.8 2.51 E-03 1.06 E—Ol

SGM 6 1 1.90 13-04 775 4.1 73.3 2.50 15-03 2.18 13-01

SGM 10 1 5.5 E-04 762 4.1 121.9 2.50 E-03 1.24 E-01

SGM 15 1 7.5 E-04 822 4.5 200.7 2.50 E-03 1.15 E-Ol

SGM 20 1 7.5 E-04 975 5.5 317.8 2.50 E-03 1.03 E-01
SFSGM 6 1 4.0 E-05 795 4.2 74.2 2.50 E-03 2.59 E-Ol
SFSGM 10 1 2.0 113-05 953 5.1 150.8 2.50 E-03 1.49 E-01

Triangular Flux (0' = 0.0018°C) - 21 Spatial Components

GM all all 1.70 E-03 9 - 54.7 1.84 13-03 8.50 E-OZ

SGM 6 1 2.20 E-03 760 4.0 137.9 1.82 E-03 1.33 E-Ol

SGM 15 1 1.25 E-03 718 4.0 328.2 1.82 E-03 8.44 E-02
SFSGM 6 1 5.50 E-05 786 4.1 141.6 1.82 E-03 1.65 E-01
SFSGM 15 1 9.50 13-05 693 3.8 314.4 1.88 E-03 9.29 E-02

 

 

163

Analysis of the triangular and step test cases with measurement errors are presented in

Table 6-2. Numerical conditions similar to the exact cases were used, a dimensionless time

step based on the sensor depth is At: = 0.06 with twenty-one nodes in the y-direction

and eleven in the x-direction for the numerical solution. Data was assumed at two-hundred
time steps and eleven components of the heat flux are estimated; twenty-two hundred heat
flux components are estimated.

The step heat flux case was more difficult to estimate. Reproducing a step change spa-
tially is considerably harder than the more gradual changing triangular case. Figures 6-6a
and b show the whole domain and SGM for r=15 estimated heat flux for the step case.
Estimates have more variability than the triangular test case. Tabulated results in the mid-
dle of Table 6-2 show that trends similar to the triangular case. The mean-squared error is
as large as 25% of the maximum heat flux for small r-values. For the triangular case the
mean-squared error did not exceed 16% of the maximum heat flux.

The whole domain (GM) approach was superior to the sequential approach (SGM and
SFSGM) for both two-dimensional test cases for r less than 10. This superior performance
is supported by the smaller mean—squared error, final column of Table 6-2, and the shorter
computational time, column seven. The mean-squared error was 7 to 83% greater for
sequential approach, compared to the whole domain approach for the triangular heat flux,
depending on the choice of r. Results were worse for the step heat flux, which had a mean-
squared error that was 17 to 144% (depending on r) greater for the sequential approach.

Computational times were adversely affected by the sequential approach. The sequential

164

a) GM
Whole Domain

  

f. ”if
\ 7 /
Z4

7

 

 

1.5\

1\
0.5\
04

Ammo—coicoEEv 55 :33

Figure 6-6a Estimated surface heat flux for step test case using data corrupted with

measurement errors (0' = 0.0025°C ). Whole domain solution

165

              
   

r=15

b) SGM

               

 

0.4

    
            

  

0.2

 

     

——

  

22>
22K”—
22 22222
2 2222222225
2%%/Z§I
£7.22: 7/j\\__
VVNK _
2: ii:
2.2/ 2 2
(V/xmyfiifl

2%
2,
22222 i
..2 2
.225255
2222 .2222... .EE'

2/

:22
. Ti/ //4»Ml...!

2/

 

15.
1
5
0.
-052
15

93222022622223 23c 28E

Figure 6-6b Estimated surface heat flux for step test case using data corrupted with

measurement errors (0' = 0.0025°C ). Sequential solution r=15

166

approach required 60 to 400% more computational time, compared to a whole domain
solution. The computational time increases proportional to the increase in the number of
future time steps, r.

The poor performance of the sequential approach for these cases can be attributed to
two related factors. The first factor is that the sequential implementation is more ill-posed
than the whole domain solution. Since the sequential solution solves over a shorter time
domain, the effect of the measurement errors is more significant and results in a more ill-
conditioned problem than the whole domain problem. Examining Table 6-2 shows that,
even though the sequential problem is more ill-posed, the magnitude of the Tikhonov
parameter in column four is decreased (except at large r—vales). Decreasing the Tikhonov

parameter is contradictory to stabilize a more ill-conditioned problem. However to obtain

the required magnitude in the sum-of-squares function, S y, the Tikhonov parameter must

be decreased. The reason for this outcome (decreasing OLT for SGM) is the second factor

in the poor performance of the sequential implementation. It is again the difficulty with the
values near the end of the sequential interval influencing the early values. Much the same
as was observed for the one-dimensional problem. Because the values near the end of the
time region are biased, the regularization parameter must be decreased to provide the
inverse solution “flexibility” to obtain the proper magnitude of the sum-of—squares func-
tion. Although decreasing the Tikhonov parameter reduces the temperature sum-of-
squares, it results in an increase in the variability of the estimated heat flux.

In contrast to the one-dimensional solution, maintaining the heat flux constant over
the sequential analysis interval does not improve the two-dimensional results. In the two-

dimensional solution it does not work because the additional stability introduced by

167

maintaining the heat flux constant over the sequential interval requires a reduction in the

Tikhonov parameter to obtain the desired sum-of-squares. The reduction in (X7- in turn

increases the variability in the estimated heat flux. If r is too large when using function
specification, the sum—of—squares can not be reduced to the desired level, regardless of the
Tikhonov parameter. The estimated results are actually worse when maintaining the heat
flux constant (with time) over the sequential interval for the two-dimensional problem.

The sequential method performs poorly for small magnitudes of r, which are in the
range that would typically be used for function specification. The estimated results for
both test cases are improved as the number of future time steps is increased. Figure 6-7
demonstrates that the effects at the end of the interval are reduced by lengthening the
sequential interval. At r=15 and 20 the mean-squared error in Table 6-2 is +1% and —12%
of the mean-squared error of the whole domain estimation, respectively. The improved
accuracy in the estimated heat flux using larger r-values is at the cost of a significantly
increased computational requirement. Sequential solutions (SGM) for larger magnitudes
of r require SOC-800% (depending on r) more computational time than the whole domain
solution. Computational time increases proportional to the increase in r. However, the
potential benefits of an on-line (real time) method may outweigh the increased computa-
tional time.

For the one-dimensional problem, and the two-dimensional problem with exact data,
an improvement in the computational speed, with a comparable accuracy in the estimated
function, was obtained when the heat flux was maintained constant over the sequential
interval, the method was denoted SFSGM. In the two-dimensional analysis with measure-

ment errors the results are not improved when the heat flux is maintained constant over the

 

 

169

sequential interval. Although computational times are comparable (to the SGM) for a
specified r—value, the error in the estimated heat flux actually increases, as seen in the final
column of Table 6-2, when maintaining the heat flux constant over the sequential interval.

It was shown with exact data that a greater number of spatial components can be esti-
mated than there are sensors available. The triangular case is reanalyzed for data with
measurement errors and results are listed in the final block of Table 6-2. The outcome was
similar to that seen in the analysis with exact data. A comparable heat flux, to that esti-
mated with few spatial components, is obtained for a whole domain solution. Using the
sequential method has a larger variability in the estimated heat flux for small r—values. The
sequential method is improved by increasing r. Computational times are significantly
greater using the sequential method.

Several approaches are possible to improve the sequential method. Before outlining
these approaches, a summary explaining the effects that are resulting in the sequential
method’s poor performance is given. It is understood that the sequential method is more
ill-posed. However, when solving with a sequential method the Tikhonov parameter is typ-
ically smaller than the whole domain solution; as r is increased the Tikhonov parameter
increases for the sequential method. The reason that the Tikhonov parameter decreases
(assuming that the residual principle is applicable), even though the sequential problem is
more ill-posed, is due to the components near the end of the time interval. These compo-
nents are difficult to estimate; when using a gradient method with the adjoint equation
approach, these components are not estimated. Instead, the adjoint function, which is pro-
portional to the computed update to the heat flux, is defined to be zero at the end of the

time interval. Consequently, the heat flux estimated at the end of the time interval is “near”

170

the initial estimate specified for the heat flux. Its deviation from the initial estimate
depends on the magnitude of the Tikhonov parameter.

As was demonstrated, lengthening the sequential interval improves the mean—squared
error in the estimated heat flux; see Table 6-2 and Figures 6-7. Lengthening the sequential
interval is also shown to increase the computational time in Table 6-2. An advantage of the
SGM, compared to the SFSGM is that functional constraints are not imposed on the solu-
tion. Consequently more than one component may be retained for a sequential interval and
the computational time can be reduced by having fewer sequential intervals. For example,
when retaining 10 components for the triangular heat flux with r=20, the computational
time is reduced to 30.8 seconds, which is about only 3% greater than the whole domain
computational time. The mean-squared error for this case, shown in the final column of
Table 6-3, is within 1% of the whole domain error; the mean-squared has increased 14%
compared to retaining only one component for r=20. Although lengthening the sequential
interval and retaining more components on time improves the computational aspects, the
process may reduce the likelihood that non-linear problems can be linearized. Since larger
sequential intervals are considered, it is less likely that the changes in thermal properties,
due to the temperature variation during the sequential interval, are negligible. This is an
important advantage of the sequential method. The validity of this assumption will depend
on magnitude of non-linearity in the problem.

The one-dimensional problem was quite insensitive to the specified initial condition.
Regardless of how the initial condition was specified, the solution required approximately
the same number of iterations to converge. This same situation was found for the two-

dimensional problem. Estimated heat flux and computational time were insensitive to the

171

procedure for setting of the initial condition for the two-dimensional problem. Whether
the initial estimates for the next sequential interval were specified from the previous inter-
val as outline in Section 5-4.2, or specified as zero, it did not significantly change the num-
ber of iterations.

A characteristic noticed for the one-dimensional problem, and two-dimensional test
cases with measurement errors, was that the Tikhonov parameter decreases for the sequen-
tial problem, compared to its magnitude for the whole domain solution. The decrease in
the Tikhonov parameter was more pronounced at smaller r-values and approached the
magnitude of the whole domain solution as r was increased. A smaller Tikhonov parame-
ter was not expected for the sequential solution because the problem is more ill-posed. It is
the sequential implementation that forces the Tikhonov parameter to be reduced, at least
for small r—values. Since the sequential solution is on the shorter time interval, the zeroth
order Tikhonov regularization is more influential because the sensitivity (or in this case
adjoint function) is not as large, particularly for r-values beyond four to five time steps
from the end of the interval, which represents a dimensionless time of 0.3. The nature of
zeroth order regularization is to bias or “drive” the estimates to zero, to reduce their vari-
ability. Consequently, to obtain the desired sum-of—squares, the Tikhonov parameter must
be decreased so as not to “drive” the estimates towards zero. As r is increased to be farther
away from this region the sensitivity (or adjoint function) is larger and the regularization
has less affect on the solution.

When solving in a sequential manner, however, there is additional information avail-
able concerning the heat flux. The heat flux need not be “driven” to zero with the zeroth

order regularization, instead it can be confined using the concept of prior information. The

172

prior information for the heat flux is specified from converged values at the previous
sequential interval. The initial estimate for the sequential interval, which is specified from

the converged estimates at the previous sequential interval, are set identically to be the

prior information. The prior information is seen to enter the solution as q pn.(r, t) in equa-

tion (4-2.4), which is the Tikhonov regularization term in the sum-of-squares function.
Hence, instead of penalizing estimates that are different from zero, it penalizes estimates
that are different from the previous estimates. For components that are far enough away
from the end of the sequential interval the prior information has little influence. Near the
end of the sequential interval, however, the prior information is more influential and helps
control components in this region.

The results using prior information are given in Table 6-3. As shown in the table, by
using the initial estimate as prior information, the sequential method is improved com-
pared to no prior information. Figure 6-8 shows the estimated triangular heat flux for r=6.
Compared to the estimated flux without prior information, see Figure 6-5b, results are sig-
nificantly smoother. For r=6 the mean-square error is reduced by 15%, and the computa-
tional time is reduced in half with prior information. For larger magnitude of r the mean-
squared error is comparable, while computational time is approximately one-half its
former value without prior information. Similar trends are shown for the step heat flux
case in Table 6-3. Notice that the magnitude of the Tikhonov parameter is greater when
using prior information, but because it penalizes changes from the prior information (esti-
mates at previous sequential interval) it does not bias the estimates as significantly. Prior
information allows for smaller r—values to be considered and reduces computational time

by requiring fewer iterations.

173

   
           

Heat flux (W/mz)

 
    

  

 

b) SGM
1i /

O.8\ , / /,/4///2 ‘

0'6 ///’/////"//////%

0\ in"? ///i///////////

$1.3) :‘Ifl, '{itgig’g/{éEg/éflé

 

Ia , I
’ ’i/
I I

 

Figure 6-8 Estimated surface heat flux for triangular test case using data corrupted with
measurement errors (6 = 0.0018°C ). Sequential solution with r=6 using prior
information.

174

Table 6-3 Estimation results for two-dimensional IHCP with simulated data corrupted
with random errors using prior information

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Analysis
domain Iterations A
Tikhonov Cpmp S y 0'5,
per time o 2
Method r I “T total seq int (sec) ( C) W/m
Triangular Flux (O’ = 0.0018°C) - 11 Spatial Components
GM all all 1.70 E—03 9 - 31.8 1.82 E-03 8.57 E-02
SGM1 6 l 1.40 E-03 396 2.1 37.6 1.80 E-03 1.16 E-01
SGMl 8 1 3.50 13-03 391 2.0 50.1 1.80 13-03 1.02 13-01
SGMl 20 10 3.00 E-03 59 3.2 19.4 1.82 E-03 8.21 E—02
SGM 20 10 1.20 E-03 90 5 30.8 1.80 E-03 8.51 E-02
SFSGMl 6 l 8.30 E-04 390 2.0 36.5 1.81 E-03 1.33 E-Ol
SFSGMI 8 1 2.50 E-03 384 2.0 48.9 1.82 E-03 1.11 E-01
Step Flux(O' = 0.0025 ° C) - 11 Spatial Components
GM all all 6.00 E-04 13 — 45.8 2.51 E-03 1.06 E-01
SGMl 6 1 9.00 E04 508 2.7 48.2 2.50 E-03 1.79 E-Ol
SGMl 10 l 3.50 E-03 423 2.3 67.8 2.51 E-03 1.30 E-01
SFSGMl 6 l 4.00 E-04 506 2.7 47.2 2.50 E-03 1.97 E-Ol
SFSGM' 10 1 3.00 E-03 388 2.1 61.8 2.52 13-03 1.19 13-01

 

 

 

 

 

 

 

 

 

 

1 Prior information used for solution

A final procedure to improve the sequential method is discussed. Although the proce-
dure is not numerically studied, it has promise for further improving the sequential
method. For all results presented herein, the Tikhonov parameter was specified to be con-

stant for the entire time region. Using the residual principle, the correct magnitude of the

175

Tikhonov parameter was identified. This procedure allowed for a direct comparison of
computational time for the whole domain and sequential methods. An alternative method
could vary the Tikhonov parameter during the sequential analysis. In regions that the heat
flux is rapidly varying, the Tikhonov parameter could be reduced to minimize the effect of
bias. Similarly, it could be increased in regions where the heat flux is relatively constant.
In a sequential implementation information from previous sequential intervals is available
to modify the magnitude of the Tikhonov parameter, as well as the magnitude of r. The

whole domain solution does not have the necessary information to select these parameters.

6-4.0 Results - Experimental Measurements

The experimental configuration previously used to estimate the (orthotropic) thermal
properties of the carbon-carbon composite is used to study the two-dimensional IHCP.
Recall that transient temperatures and heat flux were measured to estimate the thermal
properties in Chapter 3. Now the same transient temperature measurements, with the pre-
viously estimated thermal properties, are employed to estimate the two-dimensional tran-
sient surface heat flux. Surface heat flux is known for this experiment because the power to
the heater is measured. Furthermore, the experiment is symmetrically designed and heat
losses are negligible compared to the applied heat flux. Since the heater design has three
independently controlled heaters, by activating only one of the three heaters a two-dimen-
sional surface heat flux is experimentally prescribed. Additional details on the experimen-
tal configuration are discussed in Section 2-2.0 and 3-2.0.

An important issue in the study of the multi-dimensional IHCP is defining the spatial
representation of the unknown heat flux. In the analysis of simulated data, it was shown

that the gradient method could estimated more spatial components than the number of

176

sensors available. A rule used by the author is that the number of components estimated
should not exceed the number of sensors. If more components than this are estimated,
additional regularization is typically required to stabilize the solution. In the analysis of
this experimental case, no assumptions on the spatial form of the heat flux are used. It is
assumed that the heat flux is unknown only on the surface where the temperature sensors
are located. All other boundary surfaces are insulated. See Figure 3-1. Hence, the heat flux
that is estimated has a step change spatially and also a step change on time; a difficult case.
The measured surface heat flux is shown if Figure 6-9a. Note that the surface with
unknown heat flux, which is 7.62 cm long, has five sensors located above the heater at
locations, x = 0.89, 1.91, 3.18, 4.45, and 6.73 cm.

The numerical solution for this problem uses twenty-one nodes along the x-direction
for all materials. A schematic of the geometry is given in Figure 3-2. In the y-direction
there are four nodes across the mica heater, eleven nodes across the carbon-carbon, and
eleven nodes across the insulation. A numerical time step of 0.64 seconds, which is the
same as the measurement time step, is used. Because effective properties have been esti-
mated for the mica heater, the interface between the heater and specimen is modeled with
perfect contact. Although the mica heater is quite thin (0.44 mm), because its properties
include the contact resistance between the heater and the carbon-carbon, it is thermally
significant. Based on the effective properties of the mica heater, the dimensionless time
step (Fourier number) for the temperature sensors is 0.23. A dimensionless time step in
this range does not represent a difficult one-dimensional case; however, in two-dimensions

with a limited number of sensors, the problem is more difficult.

177

6000 \

5000 \

4000 \

3000 \1

2000 \

Heat flux (W/mz)

1000\

O\

 

  

 

 

 

 

80

    
 
   
 

60

t( seconds) 40

Figure 6-9a Measured surface heat for experimental case to estimate thermal
properties of carbon-carbon

178

Selecting the appropriate Tikhonov parameter requires knowledge about the errors to
apply the residual principle (Alifanov, 1994). For this two-dimensional experiment the
temperature was measured at seven locations, but two measurements are available at all
seven locations. The sensors located at the back surface of the carbon-carbon composite
provide little information compared to the sensors near the surface. Therefore the five sen-
sors located closest to the heater are used in the analysis. By comparing the two measured
signals at each location an indication of the noise in the measurements is gained. The stan-
dard deviation of the redundant measurements, from the average, is computed as

1

J, 2
62(0) = {171—122 [Yk(t)—Yk,l(t)]2} (6-4.l)

(=1

where Y_k(t) is the average temperature (for J k sensors) for sensor location k, Yk’ ,(t) is
the measured temperature for sensor location k and sensor number I , and J k is the num—

ber of sensors at location k. For this experiment J k = 2 for all k. Since so few sensors

are available, by substituting for the average temperature

Yk’1(t)+ Yk,2(t)

 

Y_t = 6-4.2
k( ) 2 ( )
equation (6-4. 1) simplifies to
l
l 2 2
67*“) = {ifl’k’ 1“)" Yk,2(t)] } (6-4-3)

Notice that this error is a function of time and sensor location k. An approximation of the
total error is obtained by averaging the error in equation (6-4.3) over time and sensor loca—

tion. The average magnitude is 0.1°C . This is about 20% larger than the measurement

179

error estimated for the one-dimensional experiment. The one-dimensional experiment had

a larger sample set, J k = 7.

The estimated heat flux is shown in Figure 6-9b and c for the whole domain and
sequential method. The two methods are remarkable similar. Table 6-4 shows results from
the analysis. The mean—squared errors are within 5% of the GM for all cases considered,
except when prior information is used. For r=6, however, the computational time 77%
greater for the SGM, which increases to 203% greater for r=10. The computational time is
reduced to the level of whole domain solution by retaining more than one component.
With r=6 and retaining two components, the computational time is cut in half and the
mean-squared error is reduced by 3%, compared to the estimates while retaining only one
component. Similar improvements are seen when retaining five components for r=10.

For a comparison, another program for estimating the surface heat flux is used to ana-
lyze the same data. The program QUENCHZD (Osman and Beck, 1989) computes surface
heat flux from internal measurements using a combined function specification regulariza-
tion method (CFSRM). In the CFSRM algorithm function specification on time and
Tikhonov regularization on space is applied. The finite element direct solver TOPAZZD is
used in this program. Heat flux estimated with QUENCHZD is shown in Figure 6-9d for
r=6. The Tikhonov parameter was specified to have the same magnitude as the whole
domain solution. Similar numerical aspects as used in the gradient method’s FCV solution
were defined for the finite element solution. Approximately the same number of nodes and
time step were specified for the numerical solution.

The estimated heat flux is again quite similar. For the CFSRM only seven spatial com-

ponents are estimated because of program limitations. Because QUENCHZD required two

180

                
   
            

     

b) GM
Whole Domain

/ , .
/ 2 2
22222222222 2
2227:2222; 2722'
222,22 ’2.— 22,//22///2//// I
2222,, 2:222 ,
2.22222. 22 .2221
22%/é l/!

2 297/why, 22'
722/2»; ,2-

 

      

      

 

\1/
S
d
n
O
C
e
S
(
t
\\\\\\1W\
00000000
mmmmmw m
6543211 1...

@533 5e 28:

Figure 6-9b Estimated surface heat flux for experimental case. Whole domain solution.

181

c) SGM
r = 10

 

 

 

 

 

 

6000 \

5000 \

000 \
2000 \

2.2.53: 522. 28:

t( seconds) 40

 

Figure 6-9c Estimated surface heat flux for experimental case. Sequential solution with

10.

r:

p—A

82

d) CFSRM
6000\ ' = 6
500m
4000.

3000\

2000‘

1000\

Heat flux (W/mz)

 

 

 

 

    

t( seconds)

Figure 6-9d Estimated surface heat flux for experimental case. Combined function
specification regularization method with r=6.

183

Table 6-4 Estimation results for two-dimensional IHCP with experimentally measured

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

heat flux
Analysis
domain Iterations A
Tikhonov Cpmp S y Gq
per time
Method r 1 “T total seq in: (sec) (°C) W/m2
Experimental Case 1020&65 - 21 Components
GM all all 1.5 E-06 l3 - 89.5 1.00 BO] 648
SGM 6 1 3.0 E-06 763 4.3 159 1.04 BO] 660
SGM 6 2 2.4 E-06 429 4.9 89.3 1.03 E-01 637
SGMl 6 1 1.3 E-OS 486 2.8 101 1.01 E-Ol 878
SGM 10 1 3.6 E-06 825 4.8 298 1.01 ED] 679
SGM 10 5 2.9 E-06 215 6.3 78.3 1.07 E—Ol 635
Experimental Case 1020&65 - 7 Components
SGM 6 l 2.6 E-08 504 3.3 104 1.04 BO] 500
SGM] 6 1 2.0 13-07 371 2.2 76.8 1.00 E-Ol 389
CFSRM 6 l 1.5 E-06 - - 7200 7.88 E-02 399
CFSRM 10 l 1.5 E-06 - - 7200 1.65 ED] 518

 

 

 

 

 

 

 

 

 

 

 

1 Prior information used

hours to obtain a solution the appropriate ocT, based on the residual principle, was not

identified. It appears the appropriate (r, (1T) for CFSRM is between the two cases listed in

Table 6-4. The estimated heat flux is not sensitive to these parameters. The biggest distinc-
tion between the two programs is the computational time. The gradient method requires on
the order of 100 seconds, while CFSRM requires on the order of 7000 seconds. These

computational times can not be directly compared because different direct solvers are

184
used. The gradient methods employ a ADI method, which typically is much faster than the

finite element method used in QUENCHZD. To understand the computational difference,
the time required to obtain one direct solution with comparable numerical conditions was
conducted. The ADI solution required 3.6 seconds while the FEM required 110 seconds; a
factor of 30. Indicating that the gradient methods are more computational efficient because
the inverse solution required a factor of 70 more computational time. Furthermore, the
CFSRM only estimated seven components, compared to 21 for the GM. As more spatial
components are estimated the CFSRM requires solving additional problems, some compu-
tational savings are carried from previous solutions, however. In contrast, the gradient
methods do not require solving additional problems when more components are consid-
ered. A very good feature for the multi-dimensional problems. It appears that the gradient
methods are more efficient for higher dimensional problems, but further cases, with simi-
lar direct solvers, or more detailed studies are need. Comparing the gradient methods with
standard function specification was not a thrust of this dissertation. Early indications are

that the gradient methods are very promising for multi-dimensional problems.

Chapter 7

SUMMARY AND CONCLUSIONS

A study of multi-dimensional inverse thermal problems was presented. Two main
problems were studied. Parameter estimation techniques were applied to characterize ther-
mal properties of a carbon-carbon composite material. The second problem developed the
sequential implementation of a gradient method, utilizing an adjoint equation approach,
for estimating the surface heat flux from internal temperature measurements. Both prob-
lems addressed multi-dimensional cases.

Thermal properties, represented by two components of thermal conductivity and the
volumetric heat capacity, were estimated for the carbon-carbon composite. One- and two-

dimensional transient experiments were analyzed to estimate the thermal properties up to

500°C . Several conclusions are supported:

1. The thermal properties of the carbon-carbon are described by a parabolic temperature

:C , and by a lin-

dependence, given in equation (3-40) and (3—41) for the k; CC and (pC)

ear temperature function, given in equation (3-37) for k: cc.
2. Experimental uncertainty in the estimated thermal properties was a maximum of 5.8%

in (pC); and 2.4% in k1“. for one-dimensional experiments. Two-dimensional

experimental uncertainty was a maximum of 6.0% in (pC); and k: CC and 4.4% in

k8

X,CC'

185

186

3. One and two-dimensional results for (p C )2 c and k; CC correlate within the experimen-
tal uncertainty.

4. An orthotropic model adequately describes the carbon-carbon material studied.

A sequential gradient method was developed to solve the IHCP. Several one and two
dimensional cases, using simulated and experimental measurements, were investigated to
characterize the sequential method. A comparison between the sequential method and a
whole domain solution was made. The following conclusions are drawn concerning the

sequential gradient method, and variations of the method that were proposed:

1. A sequential gradient method has an accuracy comparable to the whole domain gradi-
ent method, but in certain cases requires significantly more computational time.

2. In one-dimensional cases the minimum length of the sequential interval was a dimen-
sionless time of 0.3. The dimensionless time is based on the depth of the sensor nearest
the location of the estimated surface heat flux.

3. Applying function specification on time, in conjunction with the sequential gradient
method, reduced the computational requirements and allowed for sequential intervals to
be shorter than 0.3 for the one-dimensional cases. Computational requirements were
comparable for the sequential function specification gradient method and whole
domain method.

4. The two-dimensional cases required the dimensionless time to be greater that 1.0 to
obtain comparable accuracy between the sequential and whole domain solutions. The
dimensionless time is based on the depth of the sensor nearest the location of the esti-
mated surface heat flux.

5. Applying function specification on time, in conjunction with the sequential gradient
method, did not improve the computational requirements when errors were present in
the data for the two-dimensional case.

6. Retaining more than one component for a sequential interval improved computational
requirements without affecting the accuracy of the sequential gradient method for the
two-dimensional case.

7. Using prior information with the sequential gradient method for the two-dimensional
case was shown to improve the accuracy, permitted smaller sequential intervals, and
made the computational requirements competitive with the whole domain solution.

 

Chapter 8

RECOMMENDED FUTURE WORK

In the course of this dissertation several topics for future work were realized. The top-
ics are separated according to the parameter estimation or inverse heat conduction prob-
lem.

Experiments conducted to estimate the thermal properties of the carbon-carbon com-
posite material measured transient temperature and heat flux histories for discrete experi-
ments spanning the desired temperature range. Experiments with one- and two-
dimensional heat flow were analyzed independently, assuming the thermal properties were
constant for the analysis of each experiment. In the end, the discrete experiments were
combined to determine a temperature dependence for the thermal properties. One and two-

dimensional results were also compared. Topics of future direction are noted:

1. Analysis of currently available data to directly estimate temperature-dependent thermal
properties is of interest. A sequential analysis using prior information is used to accom-
plish this.

2. The optimal experiment for estimating two components of thermal conductivity and
volumetric heat capacity is needed for the case of a relatively high thermal conductivity
material (which this particular carbon-carbon has). The optimal conditions for the
design of a series of experiments is an important concept that needs attention. The
experimental series may included one- and two-dimensional experiments. As the com-
plexity of the experiments increase, more frequently several experiments may be
required to fully characterize complex materials or compound structures.

3. The model of the carbon-carbon composite can be extended to include the silicon-car-
bide layer. This case of a thin layer on a thick substrate represents an important prob-
lem. Work to estimate the thermal properties of the thin layer, in conjunction with the
underlying material, is of interest.

187

188

. Methods, both experimental and analytical, should be generalized to the three-dimen-
sional case. Developing the analysis is straight forward. Experimental techniques,
including the experimental design, are extremely important for higher dimensions and
will require significant effort.

The sequential gradient method has several areas for potential growth of the method.

. The procedure to linearized nonlinear problems in the sequential gradient solution sug-
gests large savings in computational time without significantly affecting accuracy, the
procedure needs investigating.

. Implementing on-line logic for selecting the regularization parameter, number of future
time steps, and/or the number of retained components would enhance the sequential
method. Quite frequently, the selection of these parameters is not trivial and requires
experience and an understanding of inverse problems. Since the accuracy of the solu-
tion may be sensitive to correctly selecting these parameters, the method would be
enhanced to logically select the parameter.

. In the sequential formulation presented, zeroth order Tikhonov regularization was used.
Extending the method to include first order Tikhonov regularization is appropriate.

. Prior information to spatially filter the heat flux could provide spatial regularization.

5. Implementing the methods for the three-dimensional case is of interest.

APPENDIX

APPENDIX A

FINITE CONTROL VOLUME METHOD

A-1.0 Introduction

Many methods exist for discretizing partial differential equations. Two popular meth-
ods are finite differences and finite element. An alternate approach applies conservation of
energy directly to a control volume. That method is applied to discretize the heat conduc-
tion equation here.

The problem considered is the heat conduction equation with constant thermal prop-

erties, a spatially time dependent volume energy source

 

8 er a 87' _ er a
5:0“ 5;)+$(kt5;)+g(x, y, t) — pC[§; +u0§£| (A-1.l)
8T _ ._
-k é—n—i — qi(si, t) (r --l, 2, 3, 4) (A-1.2a)
F,
OR Tlr = Ti(51" t) (i=1, 2, 3,4) (A-1.2b)
T|t=0 = T0(x, y) (A-1.2c)

The thermal conductivity is assumed orthotropic. First or second kind boundary condi—

tions are allowed on the four bounding surfaces (1“,) . Conditions may have temporal and

spatial variation. The outward normal of surface 1' is ni. Since a regular geometry is used

this is the unit normal in the x or y direction.

189

190

The FCV procedure uses energy conservation principles to derive a discretized set of
equations to solve for the temperatures. In contrast, a finite difference scheme would dis-
cretize the problem presented in equation (A-l.1) and (A-l.2), i.e. use a finite difference
approximation for the derivatives, to obtain a numerical solution. The FCV procedure is
preferred over a finite difference approach because it is more general and applies energy
conservation principles to derive the equations.

To implement the FCV procedure the spatial domain is discretized with uniform mesh

points (nodes)

(x, y) = (iAx, jAy) i: (o, 1, 2 ..... ,1), j = (o, 1, 2 ..... ,J) (11-1.3)

where the mesh lengths are

L L,
( = .5 (Ay = 72) (A-1.4)

Control volumes are the region (AxAy) surrounding the nodes with the node centrally
located. The control surfaces are at the mid-plane between the nodes. Figure A-l shows a
schematic of the nodes, control volumes, and control surfaces.

The integral energy equation (Beck et al., 1992)

J(_q . fi)dA + Jgdv = Jnggdv+ Iqu’ - fi)dA (A-1.5)
C.V C.S.

C.S. C.V.

is applied to each control volume on the domain to obtain the FCV equations. The integral
energy equation can be used to derive the describing differential equation in equation (A-
H) as well. The energy balance for a typical node in the interior is shown in Figure A-2.

Using the simplified approximation of equation (A-1.5), which assumes the quantities are

191

Control
_>| Ax I‘— Surface “Ode

 

l
_I___L__L__|____l__l___l __I____I_.__I_
| I I I | I | l I I
I O I O I O I O I O I O I O I O I O I
-—|———I-————I——-—I—-—-I---I--—-|———|—-—-I-——I— Ay
i1.1.1 I I('5J)I 11.1.1
I | l I I | | l | I
‘l"—‘T-—1"”|_—'T“"T'__|-—"I—“T—"'I— T
1 I o I o I o I o I o I o I o I o I o I
I I I I I I | I l I
'"I—__I-_T—_I__—I__T__|_—_I—_T_-I-
j=O I. 1
i=0 1 2 ..... I
Figure A-l Computational nodes and control surfaces
(1.131)
‘11.
u I'__f--‘I u
(i_1,j)o I o —I-—> o(i+l,j)
‘11: I (131’) I 41-.
.__I___.
‘11-
o
(inf—'1)

Figure A-2 Typical energy balance for an interior node

192

constant over the control volume and surface, gives

3T
_ (_ 417411 + qi+Ai+ _ 411"}: + qj+Aj,) + gV = pC-é; V + pvoui+Ai, — pvoui-Ai. (A-1.6)

The heat fluxes are taken as the mean values over the control surface and are approximated

as

(Ti,j_ Ti- l,j)

 

 

 

 

‘11; = _kx Ax (A-1.7)
qr, = -kx(Ti+1X; TU) (A-1.8)
‘11. = ‘kym’flaly- TU) (A-L9)

Assuming that the velocity is positive, the flow terms are evaluated at the upwind locations
to avoid the known difficulties associated with central differences for this term (Anderson

etal.,)
u- = C T«,- (A-l.11)

“i- = C1971- (A-1.12)

l,j

The areas are (A1.) = Ai = Ay) and (Aj. = A]. = Ax) and the volume is (AxAy). Substi-

tuting equation (A-l.7) through (A-l.12) with areas and volume into the energy balance

gives

(Tij—Ti,j—l) (Ti.j+I_Ti.j)
Ay Ax+ky Ay

(Ti.j—Ti~l.j) (Ti+1,j"Ti,j)
k1 Ax A_v+kx Ax

 

 

Ay — kv Ax

 

 

8T. .
+gAxAy = pCAxAy-a—;"1+pCv0Ay(Tl-vj— T,-_ UNA-1.13)

193

Defining the following constants

k,r k
x, = —; x = ——-”—— (A-l.14)
PCUM)2 y PC(A)’)2

and simplifying equation (A-1.13) results in an equation that represents an energy balance

on a finite control volume

V V
[(26, + A—QT" 1,1- (226x + icy“ + it'xTi+ 1, j] + [4.3713 j_ 1 — 2231131.... MT“). 1]

+ iLé = g—tTi’jm-LM)

For boundary conditions that are prescribed temperature only, after approximating the

time derivative in equation (A-1.15) the equation can be written for each interior node and

numerically solved for the all discrete temperatures. However, if flux boundary conditions

are prescribed, a FCV procedure performs an energy balance on the boundaries and results
in an equations to describe the energy balance there.

The FCV equations to represent the boundary conditions apply an energy balance at

the boundary FCV. Two possible configurations on the boundary are shown in Figure A-3.

The first is along a bounding surface but not on a comer and the second is at a comer.

The energy balance for the bounding surface node is

2 3T
—(" qi+Ai+—qjjAj_+Qj+Aj++qsth5)+gV = 9C3? V+pv0u0Ai—pvou1A,-+ (A’l°l6)

The heat fluxes (q 1' , q j, q j ) and flow term (ui ) are the same as for an interior grid point

and given in equation (A-l.8) through (A-1.10) and equation (A-l.l 1). The flow term at

the boundary surface is evaluated using the average temperature,

ui = Cp(T0,j+Tl,j)/2' The areas are A}. = A}: = Ax/2, AI. = A5 = Ay and the

194

0, '+1 0
( J ) q” (0, 1)
J o
I' _ 1 q"
ui<—I Ii? "Lb:- — 11" u!”
. 4—5 —I;I> 00’!) 3132,,- I ' €1,-
qsw‘ (0,J')| | 1* (0,0)‘ “ — " * ' ~
I I . (1,1)
L _ J qsl
q}.-
(09 j_ 1) .
a) Bounding surface node b) Comer node
x=0 (x=Qy=0)

Figure A-3 Finite control volumes along the boundary of the domain

volume is V = AxAy/Z. Substituting these relations into equation (A-l .16) and simplify-

ing gives
1 V0 , V0 1 1 '
[PMEITOJ+IZMEITLJI+ [1,181-1 -2x,r,,+>.,r,j, 11

2q52vj g0,j _ a_TOvj

—pCAx pC — a:

 

(A-1.17)
where Xx and Xx are given in equation (A-12). The prescribed heat flux (“15" j is the heat

flux on surface 32 at nodej.

Using the same procedure the following equation is derived to represent the energy

balance at a corner node

195

V V
[(423, + $970, 0 + (21', _ E1)“ 0] + {—271370, 0 + 223,10, 11

2631.1“ ‘I’ 2‘49: IA)’ 81,; 3T1.)
— + —- = — (A-1.18)
pC(AyAx) pC a:

 

The equations derived for the boundary conditions, equation (A-l7) and (A-18),

apply along boundary x = 0 and at comer (x = 0, y = 0). Similar expressions can be

derived for the remaining boundary surfaces and corners.

A-2.0 Alternating Direction Implicit Method (ADI)

Discrete equations derived can be solved using an explicit or implicit method. The
explicit is more easily implemented, however, may be unstable depending on thermal
properties and the time and space discretization. In contrast, an implicit method is more
complicated to implement, but unconditionally stable. A fully implicit method has five
unknowns for each interior node (see Figure A-2 and equation (A-1.15)) and requires
solving a penta-diagonal system for each time step. A more efficient method is the ADI
(alternating direction implicit) scheme, which requires solving two tri-diagonal system of
equations at each time step.

The time domain is uniformly discretized

t=nAt n=(l,2 ..... ,N) (A-2.l)

In the ADI method, two solutions are obtained for each time step. For the first half time
step the FCV equations are solved implicit in the x-direction and explicit in the y-
direction. For the second half time step the difference equations are implicit in the y-

direction and explicit in the x-direction. Hence, the ADI method is a two step solution. An

196

intermediate solution at time (n +1/2)At is obtained after the first step. The second

results in the solution at the end of the time step.

A-2.l ADI Equations (n+1/2)

The FCV equations using an ADI scheme proceed over the first half time step. These
equations are implicit in the x-direction and explicit in the y-direction. For illustrative pur-
poses consider the FCV equations in equation (A-l.13), (A-1.15), and (A-l.16). Introduc-
ing the forward difference for the time derivative and the implicit/explicit ADI scheme for

the first half time step produces the following equations

Interior Node: (0<i<I), (0<j<J)

 

v At v At
[(lxl' 0 {MHZ—(273+ 0 1.1-1:1/2+an+1/2:|

 

       

 

 

       

 

 

Ax 1-1,j Ax 1,} x 1+1,}
n n+l/2_
gi __J' _ (Ti) 1'}
+(ATZJ_12)tyTU+7LTU+1)p——’E_ At/2 (A2.2)
Flux B.C.(i = 0) and(0 <j< J)
n+l/2 , VON n+1/2
[I’M-Ax IT01 I‘ Ax 1'1 I
All It n+l/2 n
n n n qs,,j gO,j = (T0,j _TOJ)
Flux B.C. (i=0) and (j=0)
n+I/2 , VAI n+l/2
[I 2 ”)7.“ (”r—:1. I
An A" n n+l/2 n
n q Ax+q M T —T..
+(—2A),T0‘0+2ky7"g,l)—2 “'0 ”'0 go") (00 "1) (11-2.4)

pC(AyAx) +176 = At/2

197

The first equation is for an interior node and the next two equations are for a surface and
corner node for a prescribed flux boundary condition. For the case of temperature bound-
ary conditions only, the first equation is sufficient for numerically solving the problem.
Using a similar procedure, equations for the remaining surfaces and corner can be written.
After deriving all the FCV equations, they can be rearranged to group known information
on the left hand side of the equation and unknowns on the right. The equations represent a
set of simultaneous equations, which can be written for each j, (j = 0, 1 . . . ,J), the jth

set of equations are

I

v At
(_kx—fiyf‘kh (20. +1)+v" Ax —)T’.’j“2— 2. 1711“}: D? (O<i<I) (A-2.5b)

               

VA’ n+l/2 . vAt ““2 n .
ADMITO'I +(7nx+ 2,; LJ =00(l=0) (A-2.5a)

v At voAt
(_zlx__27)17+11/2+ (2Q + l)+v0 Ax _)TZ;1/2 = 0'; (i=1) (A-2.5c)

The right hand sides are also a function of j,

fori=0

At An 2At An 2A:
1);; _-2(1—71)T3 0+?» To 1+83,opc+‘1s1.0pcAy+q‘2’OPCAx

 

(j = O) (A-2.6a)

It At An 2At

03:2r0j_,+2(1_2)roj+2r01+1+g,1pc+qsl}.pC——A—x0<j<1 (A-2.6b)

n A__I_+An12At +An 2At

 

for0<i<l
At A 2At .
= 2(1— A)T?O+1)T?I+g:0p—C+q:l-pCAy (1:0) (A-2.7a)

198

nAI

D'.’ = 2. T’.' ._,+2(1—2I)T2I+2 I’I‘III-egII— pC

, y ,I I (O < j < J) (A-2.7b)

I: At An 2A_t_

DI.'-7IyTi,J- 1-+-2(1—---IIH)T 'J+gI'JpC+ 453i ’p—-_CAy (j: J) (A-2.7c)

for (i = I)

n n At ,. 2At . 2At
-_- 2(1—AI,)T,I0+AI,T,II+g70 " "

. p—C + qsvopCAy + q-‘b OpCAx (j = 0) (A-2.8a)

I: At An 2At
=AIT,I_ 1+2(1"}‘y)Tl.yj+)‘le+1+gt.jp_c+qsi.ijAx

 

O < j < J (A-2.8b)

n A! An 2At An 2At

D, =1. yTIJ- I+2(1- A. )TI J+gl Jp—-C+q 53’ lm+qs"1p—CA_X (j=J) (A-2.8C)

A-2.2 ADI Equations (n+1)

The FCV equations for the second-half time step for the ADI scheme proceed from
the intermediate solution at n+1/2 to the n+1 time step. For this second half step the FCV
equations are implicit in the y-direction and explicit in the x-direction (the opposite of the
first-half time step). Using the same procedure as before, the forward difference approxi-

mation is used to represent the time derivative, a simultaneous set of equations are written

for i, (i: 1,2...,I),the ith set ofequations are
20‘.» +1)T?;‘—AIT:{‘ = 133”” (j = 0) (A-2.9a)
A. Tn+1 A n+1 n+1 n+1/2
— ..v i.j- l+2( +l)TI.I. —}t yth+l = 51' O<j <J (A-2.9b)
n+1 n+1 n+1/2
4. TI,_ 1+20‘y +1)TI., = E, (j=J) (A-2.9c)

The right hand side constants are a function of 1:

forj=0

+1/2 +1/2 V A‘ +1/2
£3 =(ZII-MV)-AX)T".0 +(2I AxIT’Ifl

     

199

 

   

 

 

       

At . 2At ,. 2_A_t
+83,opC “13,, op—CAy+q:2,o CID—FAX (i=0)(A-2.10a)
n+1/2 VA! n+1/2 _AVot +1/2 +1/2
+3120 IIAé-i-‘I'I'II I-p———2CAAtII O<i<I(A-.210b)
+1/2 VA’ +1/2 _OVA’ +1/2
E3 — -(2?»I+ Ax )YJI'_ 1.x0+(2(1‘ 20
I: At An 2At An 2_A__t _
+gIOp——C+q S1. I'OpCAy+q‘4 p——-CAII(i I)(A-2.10c)
forO<j<J
n+1/2_ VOA n+1/2 VOA! n+1/2
EI. _(2(1—2.I)—E 0,; +(22I— 733—)le
I: At An 2At
+g0’jp_C+ qssz—_CAx (i: 0)(A-2.lla)
n+1/2 VON n+Il/2+ _A_Vot +1/2 n+1/2
At
+gI IpC0<z<I(A211b)
voAt _OvAt
E;+l/2 = (2kx+ Ax )T7+I1/2+ IMZIIII Ax)TJI.I+I1/2
I: At An 2At ._ I
+gIij—C+qIIIIIpCAx (z — I)(A-2.HL)
forj = J
En+l/2: 21 A n+1/2+ VON n+1/2
_[ ( - )- EAITOJ Ax T1,]

200

At & 2At ..
+g0, JpC:: :2, IJpCAx+q‘3’O pZCAAy

n+1/2_ V__.) A! MI1/2+ I+1/2+ Tn+1l2

(i: 0) (A-2.12a)

+ n 3...}. A 2A!
g‘JpC+ q'33' pCA—_y

VA’ n+1/2 VA! +1/2
nyl—LJ “(2(1-xx 2x 1,]

+g",,)At+. 2At +* _2__A1
pC q‘3’pCAy 45» pCAx

O<i<I(A-.212b)

 

n+1/2
E, (=22tI +

     

(i =I)(A2.12c)

 

A-2.3 Summary of ADI Equations

The ADI scheme results in two tri-diagonal sets of equations to solve for the tempera-

ture
{A1}I.T’.'_+,”t"+{c11T’.'I’f“2+{131}I.T;’;‘I“I2 = D? (A-2.l3)
{A2}ITZ;_II+{C2}I-T?I;l +{192}I.T§I’TI‘1 = 157*”2 (A-2.14)

where the coefficients, {Al }I-, {C l } I, {Bl }I, are given in equation (A-2.5) and D? is

given in equation (A-2.6) through (A-2.8) and coefficients {A2} 1" {C 2} I, {82} I, are

given in equation (A-2.9) and E; + “2 is given in equation (A-2. 10) through (A-2.12).

LIST OF REFERENCES

LIST OF REFERENCES

Alifanov, O. M., Artyukhin, E. A., and Rumyantsev, S. V., 1996, Extreme Methods for
Solving III-Posed Problems with Applications to Inverse Heat Transfer Problems,
Begell House Inc., New York.

Alifanov, O.M., 1994, Inverse Heat Transfer Problems, Springer-Verlag, New York.

Alifanov, O. M. and Egorov, Y. V., 1985, “Algorithms and Results of Solving the Inverse
Heat Conduction Problem in a TWO-Dimensional Formulation,” Journal Engineer-
ing Physics, Vol. 48, No. 4, 489-496.

Alifanov, O. M. and Kerov, N. V., 1981, “Determination of External Thermal Load Param-
eters by Solving the Two-Dimensional Inverse Heat-Conduction Problem,” Jour-
nal Engineering Physics, Vol. 41, No. 4, 1049-1053.

Anderson, D. A., Tannehill, J. C., and Pletcher, R. H., Computational Fluid Mechanics
and Heat Transfer, Hemisphere, New York, Chapter 4, pp. 128.

Arora, J. 8., Introduction to Optimum Design, McGraw-Hill, New York.

Artyukhin, E. A., 1996, “Account for Smoothness when Estimating Temperature-Depen-
dent Thermophysical Characteristics,” Proceeding of the Second International
Conference on Inverse Problems in Engineering: Theory and Practice, eds. D.
Delaunay, K. Woodbury, and M. Raynaud, 9-14 June 1996, LeCroisic, France,
ASME Engineering Foundation.

Artyukhin, E. A., and Gedzhadze, 1. Yu., 1994, “Sequential Regularization Solution of a
Boundary Inverse Heat Conduction Problem,” 2nd Joint Russian-America Work-
shop on Inverse Problems in Engineering, August 1994, St. Petersberg, Russia.

Balageas, D.L. and Luc, A.M., 1986, “Transient Thermal Behavior of Directional Rein-
forced Composites: Applicability Limits of Homogeneous Property Model,” AIAA
Journal, Vol. 24, No. 1, p. 109-114.

Bass, B. R., 1980, “Application of the Finite Element Method to the Nonlinear Inverse
Heat Conduction Problem using Beck’s Second Method,” ASME Journal Heat
Transfer, Vol. 102, pp. 168-176.

201

202

Baumeister, J., 1987, Stable Solution of Inverse Problems, Friedr. Vieweg and Sohn,
Braunschweig.

Beck, J. V., Blackwell, B., and Haji-Sheikh, A., 1996, “Comparison of some inverse heat
conduction methods using experimental data,” International Journal of Heat and
Mass Transfer, Vol. 39, No. 17, pp. 3649-3657.

Beck, J. V., 1996, “Parameter Estimation Concepts and Modeling: Flash Diffusivity
Application,” Proceeding of the Second International Conference on Inverse Prob-
lems in Engineering: Theory and Practice, eds. D. Delaunay, K. Woodbury, and M.
Raynaud, 9-14 June 1996, LeCroisic, France, ASME Engineering Foundation.

Beck, J. V., 1993, “Comparison of the Iterative Regularization and Function Specification
Algorithms for the Inverse Heat Conduction Problem,” Inverse Problems in Engi-
neering: Theory and Practice, eds. N. Zabaras, K. Woodbury and M. Raynaud,
ASME.Engineering Foundation.

Beck, J. V., Cole, K, Haji-Sheikh, A., and Litkouhi, B., 1992, Heat Conduction using
Green 's Functions, Hemisphere, New York.

Beck, J. V., Petrie, T. W., and Courville, G. E., 1991, “Using Parameter Estimation to Ana-
lyze Building Envelope Thermal Performance,” In-Situ Heat Flux Measurements
in Buildings - Applications and Interpretation of Results, Steven N. Flanders, Edi-

tor, US Army Cold Regions Research and Engineering Laboratories, Hanover,
NH., Special Report 91-3, pp. 161-191.

Beck, J .V., and Osman, A.M., 1991, “Sequential Estimation of Temperature-Dependent
Thermal Properties,” High Temperature-High Pressure, Vol. 23, pp. 255-266.

Beck, J. V. and Murio, D. A., 1986, “Combined Function Specification-Regularization
Procedure for Solution of Inverse Heat Conduction Problem,” AIAA, Vol. 24, No.
1, Pp. 180-185.

Beck, J .V., Blackwell, B.,and St. Clair, CR, 1985, Inverse Heat Conduction, Wiley, New
York.

Beck, J .V. and Arnold, K.,l977, Parameter Estimation in Engineering and Science, Wiley,
New York.

Bines, E. B., 1993, “The Protection of Carbon-Carbon Against Oxidation at Elevated
Temperature,” Essential of Carbon-Carbon Composites, ed. C.R. Thomas, Royal
Society of Chemistry, Cambridge, p. 228.

Busby, HR. and Trujillo, D.M., 1985, “Numerical Solution to a 'I‘wo-Dimensional Inverse
Heat Conduction Problem,” International Journal for Numerical Methods in Engi-
neering, Vol. 21, pp. 349-359.

203

Chamis, CC, 1973, “Computerized Multilevel Analysis for Multilayer Fibre Compos-
ites,” Computers and Structures, Vol. 3, pp. 467—482.

Dowding, K., and Beck, J. V., 1994, “Measurements of Transient Temperatures and Heat
Fluxes in a Composite Material for the Estimation of Thermal Properties - Part 11:
Property Determination and Heat Flux Estimation,” Michigan State University,
MSU-ENGR-004-94, East Lansing, MI.

Emery, AF and Fadale, T. D., 1996, “Design of Experiments Using Uncertainty Inforrna-
tion,” ASME Journal of Heat Transfer, Vol. 118, pp. 532-538.

Fadale, T. D., Nenarokomov, A. V., and Emery, A.F., 1995a, “Two Approaches to Optimal
Sensor Locations,” ASME Journal of Heat Transfer, Vol. 117, pp. 373-379.

Fadale, T. D., Nenarokomov, A. V., and Emery, A.F., 1995b, “Uncertainties in Parameter
Estimation: the Inverse Problem,” International Journal of Heat Transfer and
Mass Transfer, Vol. 38, no. 3, pp. 511-518.

Gamier, B., Delaunay, D.and Beck, J .V., 1993, “Measurement of Surface Temperature of
Composite Materials for the Optimal Estimation of Their Thermal Properties,”
Proceedings of the 13th European Conference on Thermophysical Properties, Lis-
bon, Portugal.

Gamier, B., Delaunay, D., and Beck, J .V., 1992, “Estimation of Thermal Properties With-
out Instrumentation Inside the Samples,” International Journal of Thermophysics,
Vol. 13, No.6, pp. 1097-1111.

Guo, L. and Murio, D., 1991, “A Mollified Space-Marching Finite-Difference Algorithm
for the TWO-Dimensional Inverse Heat Conduction Problem with Slab Symmetry,”
Inverse Problems, Vol. 7, pp. 247-259.

Haji-Sheikh, A. and Buckingham, ER, 1993, “Multidimensional Inverse Heat Conduction
Using the Monte Carlo Method,” ASME Journal Heat Transfer, Vol. 115, pp. 26-
33.

Han, LS. and Cosner, A.A., 1981, “Effective Thermal Conductivities of Fibrous Compos-
ites,” ASME Journal of Heat Transfer, Vol. 103, pp. 387-392.

Harris, J.P., Yates, B., Batchelor, J ., and Garrington, P.J., 1982, “The Thermal Conductiv-
ity of Kevlar Fibre-Reinforced Composites,” Journal of Materials Science, 17, pp.
2925-2931.

Hensel, E., 1991, Inverse Theory and Applications for Engineers, Prentice-Hall, New Jer-
sey.

204

Hensel, E. and Hills, R., 1989, “Steady-State Two-Dimensional Inverse Heat Conduction,”
Numerical Heat Transfer; Part B, Vol. 15, pp. 227-240.

Hsu, T., Sun, N. G., Chen, G. G., and Gong, Z., 1992, “Finite element formulation for two-
dimensional inverse heat conduction analysis, ASME Journal Heat Transfer, Vol.
114, pp. 553-557.

Imber, M., 1975, “'I‘wo-Dimensional Inverse Conduction Problem - Further Observations,”
AIAA Journal, Vol. 13, No. 1, pp. 114-115.

Imber, M., 1974, “Temperature Extrapolation Mechanism for Two-Dimensional Heat
Flow,” AIAA Journal, Vol. 12, No. 8, pp. 1089-1093.

Incropera, F. P. and DeWitt, D. P., 1990, Introduction to Heat Transfer, Second Edition,
Wiley, NY.

Ingham, DB. and Yuan, Y., 1994, The Boundary Element Method for Solving Improperly
Posed Problems, Computational Mechanics, Southampton, UK.

J arny, Y. C. and Beck, J. V., 1995, “Tutorial Solution of a Heat Conduction Parameter Esti-
mation Problem Using the Adjoint Method,” Proceeding of Seventh Inverse Prob-
lems in Engineering Seminar, Busby, H. R. and Guo, L. eds., Ohio State
University, Columbus, OH.

J amy, Y., Ozisik, N. M., and Bardon, J. P., 1991, “A General Optimization Method Using
Adjoint Equation for Solving Mulitdimensional Inverse Heat Conduction,” Inter-
national Journal Heat and Mass Transfer, Vol. 34, No.11, pp. 2911-2919.

Kerov, N.V., 1983, “Solution of the TWO-Dimensional Inverse Heat Conduction Problem
in Cylindrical Coordinate System,” Journal of Engineering Physics, Vol. 45, No. 5,
pp. 1245-1249.

Kulkami, M. R. and Brady, R. P., 1997, “A Model of Global Thermal Conductivity in
Laminated Carbon/Carbon Composites,” Composites Science and Technology, Vol.
57, No. 3. Pp. 277-285.

Kurpisz, K. and Nowak, A.J., 1995, Inverse Thermal Problems, Computational Mechan-
ics, Southampton, UK.

Lamm, P. K., 1995, “Future-Sequential Regularization Methods for Ill-Posed Volterra
Equations: Applications to the Inverse Heat Conduction Problem,” Journal Math.
Analysis and Applications, Vol. 195, pp. 469-494.

Lamm, P. K., 1993, “Inverse Problems and Ill-posedness,” Proceeding of the First Intema-
tional Conference on Inverse Problems in Engineering: Theory and Practice, eds

205

N. Zabaras, K. Woodbury, and M. Raynaud, 13-18 June 1993, Palm Coast, FL,
ASME Engineering Foundation, pp. 1-10.

Lamm, P. K., 1990, “Regularization .and the Adjoint Method of Solving Inverse Prob-
lems,” Lectures given at 3rd Annual Inverse Problems in Engineering Seminar,
Michigan State University, East Lansing, MI, 25-26 June 1990.

Lawson, CL, 1974, Solving Least Squares Problems, Prentice-Hall, Englewood Cliffs,
New York.

Lee, H,J., and Taylor, RE, 1975, “Thermophysical Properties of Carbon/Graphite Fibers
and MOD-3 Fiber-Reinforced Graphite,” Carbon, Vol. 13, pp. 521-527.

Loh, M. and Beck, J. V., 1991, “Simultaneous Estimation of Two Thermal Conductivity
Components from Transient Two-Dimensional Experiments,” paper 91-WA/HT-l 1
presented at ASME Winter annual meeting, Atlanta, GA.

Loulou, T., Artyukhin, E. A., and Bardon, J. P., 1996, “Estimation of the Time-Dependent
Thermal Contact Resistance at the Mold-Casting Interface,” Proceeding of the
Second International Conference on Inverse Problems in Engineering: Theory and
Practice, eds. D. Delaunay, K. Woodbury, and M. Raynaud, 9-14 June 1996,
LeCroisic, France, ASME Engineering Foundation.

Luenberger, D. G., 1973, Introduction to Linear and Nonlinear Programming, Addison-
Wesley, Massachusetts.

Minkowycz, W.J., Sparrow, E.M., Schneider, GE, and Pletcher, R.H.(eds.), 1988, Hand-
book of Numerical Heat Transfer, Wiley Science, New York.

Murio, D.A., 1993a, The Mollification Method and the Numerical Solution of Ill-Posed
Problems, Wiley, New York

Murio, D., 1993b, “On the Numerical Solution of the Two-Dimensional Inverse Heat Con-
duction Problem by Discrete Mollification,” Proceeding of the First International
Conference on Inverse Problems in Engineering: Theory and Practice, eds. N.
Zabaras, K. Woodbury, and M. Raynaud, 13-18 June 1993, Palm Coast, FL,
ASME Engineering Foundation, Florida, pp. 17-21. .

Osman, A. M., Dowding, K. J ., and Beck, J. V., 1997, “Numerical Solution of the General
TWO-Dimensional Inverse Heat Conduction Problem (IHCP),” ASME Journal of
Heat Transfer, Vol. 119, pp. 38-45.

Osman, A. M. and Beck, J. V., 1990, “Investigation of Transient Heat Transfer Coeffi-
cients in Quenching Experiments,” ASME Journal Heat Transfer, Vol. 112, pp.
843-848.

206

Osman, A.M. and Beck, J.V., 1989, “Nonlinear Inverse Problem for the Estimation of
Time-and-Space Dependent Heat Transfer Coefficients,” Journal of Thermophys-
ics and Heat Transfer, Vol. 3, No. 2, pp. 146-152.

Osman, A.M. and Beck, J.V., 1989, “QUENCHZD: A General Computer Program for
TWo-Dimensional Inverse Heat Transfer Problems,” Heat Transfer Group, Depart-

ment of Mechanical Engineering, Michigan state University, Report. MSU-ENGR-
89-017, October 1, 1989.

Ozisik, N. M., 1993, Heat Conduction, 2nd edition, Wiley Science, New York.

Pasquetti, R., and Le Niliot, C., 1991, “Boundary Element Approach for Inverse Heat
Conduction Problems: Application to a Bidimensional Transient Numerical Exper-
iment,” Numerical Heat Transfer; Part B, Vol. 20, pp. 169-189.

Raynaud, M. and Beck J. V., 1988, “Methodology for Comparison of Inverse Heat Con-
duction Methods,” ASME Journal of Heat Transfer, Vol. 110, pp. 30-37

Reinhardt, H. J. and Hao, N. H., 1996a, “A Sequential Conjugate Gradient Method for the
Stable Numerical Solution to Inverse Heat Conduction Problems,” Inverse Prob-
lems, Vol. 2, pp. 263-272.

Reinhardt, H. J. and Hao, N. H., 1996b, “On the Numerical Solution of Inverse Heat Con-
duction Problems by Gradient Methods,” Proceeding of the Second International
Conference on Inverse Problems in Engineering: Theory and Practice, eds. D.
Delaunay, K. Woodbury, and M. Raynaud, 9-14 June 1996, LeCroisic, France,
ASME Engineering Foundation.

Savage, 6., 1993, Carbon-Carbon Composite, Chapman and Hall, London.

Scott, ER, and Beck, J. V., 1992a, “Estimation of Thermal Properties in Carbon/Epoxy
Matn'x Materials During Curing,” Journal of Composite Materials, Vol. 26, No. 1,
pp. 21-36.

Scott, ER, and Beck, J. V., 1992b, “Estimation of Thermal Properties in Epoxy Matrix/
Carbon Fiber Composite Materials,” Journal of Composite Materials, Vol. 26, No.
1, pp. 132-149.

Scott, E. P. and Beck, J.V., 1989, “Analysis of Order of the Sequential Regularization
Solutions of Inverse Heat Conduction Problems,” ASME Journal of Heat Transfer;
Vol. 111. Pp. 218-224.

Shapiro, A. B., 1986, “TOPAZ2D a TWO-Dimensional Finite Element Code for Heat
Transfer Analysis, Electrostatic, and Magnetostatic Problems,” Lawrence Liver-
more National Laboratory, CA.

207

Taktak, R., Beck, J .V., and Scott, E., 1993, “Optimal Experimental Design for Estimating
Thermal Properties of Composite Materials,” International Journal of Heat Mass
Transfer, Vol. 36, No. 12, pp. 2977.

Taktak, R., 1992, “Design and Validation of Optimal Experiments for Estimating Thermal
Properties of Composite Materials,” Ph.D thesis, Michigan State University, East
Lansing, Michigan.

Taylor, R.E., Jortner, J ., and Groot, H., 1985, “Thermal Diffusivity of Fiber-Reinforced

Composites using the Laser Flash Technique,” Carbon, Vol. 23, No. 2, pp. 215-
222.

Tikhonov, AN. and Arsenin, V.Y., 1977, Solutions of Ill-Posed Problems, Winston and
Sons, Washington DC.

Truffart, B., Jamy, Y., and Delaunay, D., 1993, “A General Optimization Algorithm to
Solve 2-d Boundary Inverse Heat Conduction Problems Using Finite Elements,”

Proceedings of Inverse Problems in Engineering: Theory and Practice, ASME,
Florida, pp. 53-60.

Tseng, A.A. and Zhao, F. 2., 1996, “Multidimensional Inverse Transient Heat Conduction
Problems by Direct Sensitivity Coefficient Method using a Finite-Element
Scheme,” Numerical Heat Transfer, Part B, Vol. 29, pp. 365-380.

Tseng, A.A., Chen, TC, and Zhao, P. Z., 1996, “Direct Sensitivity Coefficient Method for
Solving Two-Dimensional Inverse Heat Conduction Problems by Finite-Element
Scheme,” Numerical Heat Transfer, Part B, Vol. 27, pp. 291-307.

Tuan, P.-C., Ji, C.-C., Fong, L.-W., and Huang, W.-T., 1996, “An Input Estimation
Approach to On-Line 'I‘wo-Dimensional Inverse Heat Conduction Problems,”
Numerical Heat Transfer; Part B, Vol. 29, pp. 345-363.

Ulbrich, A., Beck, J. V., and Dowding, K., 1993, “Measurements of Transient Tempera-
tures and Heat Fluxes in a Composite Material for the Estimation of Thermal Prop-
erties - Part 1: Experimental Aspects,” Michigan State University, MSU-ENGR-
009-93, East Lansing, MI.

Zabaras, N. and Yang, G., 1996, “Inverse Design and Control of Microstructural Develop-
ment in Solidification Processes with Natural Convection,” presented at the 3lst
National Heat Transfer Conference Houston, Texas, 3-6 August, 1996, edited by V.
Prasad et al., ASME, New York.

Zabaras, N. and Liu, J., 1988, “An Analysis of Two-Dimensional Linear Inverse Heat
Transfer Problems Using an Integral Method,” Numerical Heat Transfer, Vol. 13,
pp. 527-533.

 

"{1111111111lll'll'llf