NONEQUILIBRIUM THERMODYNAMICS
0F FLUID MIXTURES:
PRINCIPLES, PERTURBATION METHODS,
AND PARAMETER ESTIMATION

Thesis for the Degree of Ph. D.
MICHIGAN STATE UNIVERSITY
JOHN LESTER BARTELT
1968

mm

F) LIBRA 1’" '

I a
§«ii.;"rxi.;;.r: Static ‘

. {scram-sir)! (
Lrantrtam-r

This is to certify that the

thesis entitled

NONEQUILIBRIUM THERMODYNAMICS OF FLUID MIXTURES:
PRINCIPLES, PERTURBATION METHODS, AND PARAMETER

ESTI MAT I ON
presented by

John Lester Bartelt

has been accepted towards fulfillment

of the requirements for

Ph.D. degree in Chemistry

Wflw

Major professor
Frederick H. Horne

 

Date ll November 1968

 

0-169

I lmomn or i
I ms a sour
: angx mm ma

Ll ”IV .IInI-E

ABSTRACT

NONEQUILIBRIUM THERMODYNAMICS OF FLUID MIXTURES:
PRINCIPLES, PERTURBATION METHODS,
AND PARAMETER ESTIMATION

BY

John Lester Bartelt

A fundamental nonequilibrium thermodynamic theory
and powerful, accurate methods of exploiting it have been
developed to deal with the intricate problems encountered
when various, interacting transport phenomena occur simul-
taneously. Critical examination of the fundamental theory
has led to a reformulation which embodies a new approach
to the specification of phenomenological relations through
the rational use of the entropy inequality and proper
enumeration of independent variables. The approach re-
solves most of the ambiguities in the earlier theory with
regard to validity of a Gibbsian differential equation in
nonequilibrium situations, pressure and mechanical equi-
librium criteria, kinetic energy of diffusion, and the
appropriateness of the Newtonian stress formula for mix-
tures. A perturbation expansion method has been developed

to facilitate the solutions to a wide range of problems in

John Lester Bartelt

transport theory. Coupled with the weighted residual method
and finite transform method, this approach is shown to be an
extremely versatile technique for handling such difficult
problems as the solution to partial differential.equations
with variable coefficients and the convection problem of
pure thermal diffusion. Some practical formulas are de-
rived for the latter case. A detailed theory of estimation
of nonlinearly interacting parameters from data with cor-
related errors is presented together with a computational
algorithm which has distinct advantages over the ordinary

least squares method.

NONEQUILIBRIUM THERMODYNAMICS OF FLUID MIXTURES:
PRINCIPLES, PERTURBATION METHODS,

AND PARAMETER ESTIMATION

BY

John Lester Bartelt

A THESIS

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

DOCTOR OF PHILOSOPHY

Department of Chemistry

1968

To Brooks

ii

ACKNOWLEDGMENTS

I would like to express my thanks to the Dow
Chemical Company, the Lubrizol Foundation, and the
National Institutes of Health for the fellowships which
provided my financial support.

I am grateful to my colleagues Dr. T. G. Anderson
and Mr. V. A. Nicely for the many hours of fruitful dis-
cussion and for the numerous suggestions and comments
which helped to make this thesis a reality. I would also
like to thank Mr. D. Knirk for making his plotting program
available and for his interest in the computer aspects of
this work.

I wish especially to express my sincere appreci-
ation to Dr. F. H. Horne, whose counsel and guidance have
been of immeasurable value. He has been encouraging,
critical, and above all a conscientious teacher.

Finally, I am deeply grateful to my wife Brooks,
who endured my difficult times, and whose understanding

and encouragement were unlimited.

iii

TABLE OF CONTENTS

ACKNOWLEDGMENTS O O O O C O O O O O C C 0

LIST OF FIGURES . . . . . . . . . . . . .

Chapter

I.

II.

III.

Iv.

V.

INTRODUCTION 0 O O I O O O O O O O

A THERMODYNAMIC THEORY OF FLUID MIXTURES

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

The Clausius-Duhem Inequality

Thermodynamic Processes . . .

Constitutive Assumptions for a

Fluid Mixture . . . . . . .
6. Consequences of the Entropy

Inequality . . . . . . . . .

7. Ordinary Binary Fluid Mixtures

(The Linear Theory) . . . .

8. The Transport Equations . . .

9. Discussion . . . . . . . . . .

£11wa
0 O

PERTURBATION EXPANSION METHODS . .

1. Introduction . . . . . . . . .
2. The Experiment . . . . . . . .
3. Pure Thermal Diffusion
Transport Equations . . . .
4. The Perturbation Expansions .
5. The Perturbation Equations . .
6. Summary . . . . . . . . . . .

PARAMETER ESTIMATION . . . . . . .

Introduction . . . . . . . .
. Generalized Least Squares
Adjustment of Data . . . .

NH

3. An Alternative Approach . . .

4. Fringe Shape Analysis . .

5. Time Constant Analysis . .

6. Summary . . . . . . . . . . .
CONCLUSION . . . . . . .

iv

Kinematics and Equations of Ba ances

Page
iii

vi

18
26
3O

31
34
38
51
56
59

59
65

69
72
75
96

97

97
100
107
111
115
117

118

BIBLIOGRAPHY

APPENDIX A.

APPENDIX B.

APPENDIX C.

RESTRICTIONS ON PHENOMENOLOGICAL
COEFFICIENTS IMPOSED BY POSITIVE
ENTROPY PRODUCTION

PARAMETER ESTIMATION PROGRAMS

VELOCITY AND TEMPERATURE PROFILES
FOR THE CONVECTIVE HEAT LOSS

PROBLEM

Page

122

125

130

148

LIST OF FIGURES
Figure Page

1. Vertical velocity as a function of r at
the center of the cell (2 = 0.5) . . . . . . 150

2. Radial velocity as a function of
r at z = 0.25. O O I O O O O C O O O O O O O 151

3. Temperature correction versus r at z = 0.5
for (Tr - Tm + 1/2 AT) = - 5.0 o o o o o o o 154

4. Temperature correction versus r at z = 0.5
for (Tr - Tm + 1/2 AT) = 0.0 o o o o o o o o 155

5. Temperature correction versus r at z = 0.5
for (Tr - Tm + 1/2 AT) = 5.0 o o o o o o o o 156

vi

CHAPTER I
INTRODUCTION

The mathematician has often been described as one
who neither knows what he is talking about nor cares
whether what he says is true. The physical scientist
often prides himself on his practical understanding,
thoroughness, and insatiable curiosity. Between such
extremes of abstraction and realism it would hardly
seem possible that there should be much progress. Yet
the trend has been quite otherwise and today the sci-
entist is increasingly aware of his need for mathematical
insight and the mathematician proves more and more the
stimulation of physical problems. The progress of the
moment, at any rate in the science of materials, lies in
the regions where mathematics and physical science have
common interests.

If the mathematician has little care whether what
he says is true, it is only in the sense that his primary
concern is with inner consistency and deductive conse-
quences of an axiomatic theory. He is content with cer-
tain undefined quantities and his satisfaction lies in
the structure into which they can be built. If the sci-

entist regards himself as dedicated to the progression

of understanding, then he cannot rest content with par-
ticular details. It is his understanding of the common
features of diverse problems that allows him to progress
and hence he must be concerned with abstraction and gen-
eralization. It is the business of mathematical theory
to provide just such an abstraction and generalization,
but it will do it in its own fashion and use the axio-
matic method. From what at first seem rather farfetched
abstractions and assumptions, it will produce a coherent
body of consequences. Insofar as these consequences cor-
respond to the observable behavior of the materials the
scientist handles, he will have confidence in the mathe-
matical theory and its foundations. The theory itself
will have been used to design the critical eXperiments
and to interpret their results. If there is complete
discordance between the valid expectations of the theory
and the results of critically performed experiments, the
theory may be rejected. Some measure of disagreement may
suggest modification of the theory. Agreement within the
limits of experimental error gives confidence in the math-
ematical model and opens the way for further progress.
In particular nonequilibrium thermodynamics provides a
model of the real world in which the scientist can have
a high degree of confidence.

Although the early steps in the development of the

theory of irreversible processes took place over a hundred

years ago, the cornerstone of the modern theory was laid
in 1931 by Onsager. Over the past twenty-five years a
coordinated theory based on the fundamental work of
Onsager has been developed by Prigogine (1947, 1955),

de Groot (1945, 1955), de Groot and Mazur (1962), Meixner
and Reik (1959), Kirkwood and Crawford (1952), Fitts
(1962), and many others. It is also clear that over the
last decade there has been a renaissance of interest in
thermodynamics and rational mechanics in the mathematical
world. It has been fairly widespread and attracted the
attention of many mathematicians whose abilities are of
the first order. If a few names are to be singled out,
it would probably not be unfair to the others to select
those of Truesdell and Toupin (1960) and Coleman and N011
(1963); their scholarship and extensive writing on the
foundations of continuum mechanics have been of great
influence.

Even though the development of the theory has been
extensive, the fact remains that, because of the diversity
of phenomena and range of material behaviors that come
within the scope of nonequilibrium thermodynamics, it is
unlikely that every possible area for study has been
examined. The investigation of one such area initiated
most of the research in this thesis: convection in pure

thermal diffusion experiments.

Convection in liquids due to temperature, density,
or composition gradients is an extremely delicate phenomenon.
The pure thermal diffusion experiments of Thomaes (1951)
and of Agar and Turner (l960),ampng others, indicate that
proper interpretation of their results can only be made if
the convection that is apparently present can be accounted
for. Because of its simplicity and sensitivity the pure
thermal diffusion method has the potential, if the convec—
tion problem can be solved or eliminated, of being an impor-
tant source of transport parameters of interest not only to
physicists and chemists but engineers and biologists as
well. Unfortunately, it appears that the rule of thumb
in the past has been to neglect convection and its effects
whenever possible. Because of the paucity of theoretical
and practical work on convection (hydrodynamicists had
never considered convection under conditions appr0priate
to the thermal diffusion experiments), it was impossible
to establish the convective behavior of a fluid mixture in
a thermal diffusion experiment without a re-examination of
the theory and without new (or better) mathematical methods.

Besides convection, there were and are many other
unsolved transport problems. Throughout the research which
has led to this thesis, the emphasis has been on the devel-
opment of a generalized theory and of generalized methods
of using it. Most of the ambiguities of the previous

theory have now been resolved, and powerful methods have

been obtained for applying the theory to particular experi-
ments. Although examples are given only for pure thermal
diffusion, the theory and methods are sufficiently general
that a great variety of other problems may be solved.

A reformulation of the nonequilibrium thermo-
dynamic theory of fluid mixtures is presented in Chapter
II. The questions of theoretical interest which are dealt
with in this chapter are: (l) validity of a Gibbsian dif-
ferential equation in nonequilibrium situations; (2)
neglect of inertial and viscous terms in the phenomeno-
logical equations; (3) neglect of the kinetic energy of
diffusion; (4) definition of mechanical equilibrium; and
(5) appropriateness of the Newtonian stress tensor for
mixtures. The resolution of each of these problems is
accomplished in this chapter through the development of
a phenomenological theory for "ordinary" fluid mixtures.

Chapter III is concerned with the development of
the methodology required to solve the extremely complex
transport equations deduced in Chapter II. In particular
it is shown that the perturbation expansion method is an
extremely versatile technique for handling such difficult
problems as the solution to partial differential equations
with variable coefficients and the convection problem of
pure thermal diffusion. Some practical formulas are pre-
sented for the interpretation of pure thermal diffusion
experiments including the possibility of convection and

variable phenomenological coefficients.

CHAPTER II

A THERMODYNAMIC THEORY OF FLUID MIXTURES

1. Introduction

The success of the thermodynamic theory of irre-
versible processes in its application to a wide variety
of problems in the physics and chemistry of materials has
led to a renewal of interest in the fundamentals of the
theory. The material presented in this chapter is in fact
a result of the need for better understanding of the prin-
cipal assumptions and approximations made in the theory.
The purpose here is not the derivation of a new theory; it
is instead a reformulation in which the nature of these
assumptions is clarified. Moreover, the reformulation has
uncovered many new features of the theory which previously
were buried in assumptions too casually accepted.

There is an extensive literature for the thermo-
dynamics of irreversible processes (TIP) or nonequilibrium
thermodynamics (NET) and its applications. Although a
literature survey is not intended here, a reference to TIP
or NET as it is currently practiced is generally meant to
be a reference to the works of Fitts (1962), de Groot and
Mazur (1962), Hasse (1963), and Meixner and Reik (1959).

TNIere is in addition a large volume of material in the

technical journals, much of which can probably be found
most easily by consulting the extensive bibliography of
Fitts; specific articles of particular importance from
the literature will be cited individually as the need
arises.

The objectives of this chapter are (l) the develop-
ment of a set of macrosc0pic transport equations from basic
principles which describe the behavior of a liquid mixture,
with an attempt at each stage of the develOpment to indi-
cate clearly the introduction of any necessary assumptions
or approximations; and (2) the comparison of the differ-
ences between the results obtained here and those of the
earlier theory (see Horne, 1966). In View of these objec-
tives it will be necessary in the remainder of this section
to give a concise exposition of TIP and to indicate at each
opportunity the need for a more rigorous or objective ap-
proach. There is no doubt that in such a short exposition
the case for TIP may not be presented as fairly as possible
and where matters of Opinion are involved it might be wise
to consult the original sources.

Typically, TIP begins with the fundamental princi-
ples of conservation of mass, momentum, and energy. These
basic postulates are usually given in integral form, and,
in regions where the field variables are sufficiently smooth,
the integral balance or conservation equations may be re-

placed by local differential equations. The conservation

of mass condition leads to the differentiad equation

usually called the equation of continuity,

api

ss— + Y ‘(9191’ = 0 ' (1'1)

where pi is the local density of component 1, Bi is the
velocity of the local center of mass of component i, and
it assumed that there are no sources of mass (i.e., chemi-
cal reactions). The conditions of conservation of linear

and angular momentum lead to an equation of motion for the

fluid,

399

55— + V - (puu) = pX + V - c (1.2)
where p is the total local density,
r
p E zpi , (1.3)
i=1

and u is the local center of mass or convective velocity

p.u. (1.4)
~ . 1 1~1

J.

'O
:2
III

II MB

The quantity X in (1.2) is the external force per unit
mass acting on the fluid system, and the tensor 0 is the
stress tensor, the negative of the pressure tensor, which

gives the surface force acting on a fluid element. Finally,

conservation of energy is expressed by

———— + y - J = 0 , (1.5)

where E is the total specific energy and J is the total

T
flux of the total energy.

E

It should be mentioned that the theory is already
informationally deficient. Unless the external force, the
stress tensor, and the total energy flux are independently
specified there are more unknowns than equations. The pro—
cedure in TIP is to provide phenomenological relations
(sometimes called constituitive relations) for these quan-
tities, although the origin of these phenomenological rela-
tions is not determined with certainty. In some formula-
tions of TIP a prescription is given for obtaining these
phenomenological equations as will be seen later, but it
is not clear that these prescriptions are more than arbi-
trary. The develOpment given in this chapter approaches
the determination of phenomenological relations in a more
systematic if not more reasonable fashion.

Another source of consternation is that TIP generally
fails to consider the possibility of determining equations
of motion for the individual components of the mixture.
Clearly momentum is not conserved for each component of the
mixture since there is exchange between components. This
is not a disadvantage however, since it will be shown that
equations of motion for the individual components can be

formulated by the use of componentwise or partial stress

10

tensors without loss of information and with great increase
in the sc0pe and the generality of the theory. This single
aspect has led to many of the new results produced in this
reformulation.

At this point TIP makes a transition from the hydro-
dynamic stage represented by equations (l.l), (1.2), and
(1.5) to the thermodynamic stage. This is usually accom-
plished by introducing the assumption of local equilibrium.
For example Fitts (1962) uses the following postulate:

For a system in which irreversible processes are
taking place, all thermodynamic functions of state
exist for each element of the system. These thermo-
dynamic quantities for the nonequilibrium system are
the same functions of the local state variables as
the corresponding equilibrium thermodynamic quantities.

The concepts defining the thermodynamic functions of state
are not generally so restrictive as to exclude the possi-
bility of extending the definitions to nonequilibrium situ-
ations. It is another matter, however, to assume that the
same functional relations exist outside of the domain for
which they were established. Some attempt should be made
to investigate the resonableness of this claim. Prigogine
(1949) in fact made such an investigation and found that
the postulate was reasonable for mixtures of monatomic
gases sufficiently close to equilibrium. Little else is
known about the validity of the postulate except by com—

parison of the results which come from it with experiment.

In this respect the postulate has performed very well for

11

systems whose gradients of the thermodynamic functions are
small and evolve slowly.

Despite these practical successes, there are still
some inconsistencies which require clarification. One in
particular was the stimulus for much of this research.
According to the postulate of local equilibrium the spe-
cific internal energy, E} is related to the specific

entropy, S, by the equation of Gibbs;

r

+ 2 u.d
i=1 1

Di

'5— r (1.6)

 

 

d§=Td§-pd(%-

 

where

 

 

is the specific volume, “i are the specific

0 0
chemical potentials, and (7% are concentrations (mass

 

fractions). With regard to the quantity p Fitts makes
the statement:
For a nonequilibrium system, the pressure p (in equa-
tion (1.6)) is the same function p(p,E) as for an
equilibrium system. In this case, p is not the pres-
sure in the usual sense; i.e., p is not the normal
force per unit are exerted by the fluid.
In hydrodynamics the pressure or force per unit area is
given in terms of the trace of the stress tensor. It
seems strange that after stating a postulate to the con-
trary, Fitts reaches the conclusion that the quantity p
in equation (1.6) is no longer the pressure normally used
while the temperature, chemical potentials, etc. all cor-

respond to the concepts typically associated with them.

Furthermore TIP in general fails to discuss the correspondence

12

between the quantity p and the trace of the stress tensor.
This makes any critical analysis of the criterion for
mechanical equilibrium within the context of TIP nearly
impossible.

At this juncture of the thermodynamic stage of TIP
an entropy balance equation is deduced from the energy
equation (1.6). It is always possible to write this bal-
ance equation in the form

So. . ¢
32"“? :Ts‘r' ”-7)

where Js is the entrepy flux and ¢ is the internal entr0py
production. It is important to note that the separation
of terms in this equation between a flux and source is
fairly arbitrary. Furthermore the introduction of this
equation to the list of transport equations does not in-
crease its information content. The real need for this
equation is for application of the second law of thermo-

dynamics. The form of the second law most appropriate

for nonequilibrium thermodynamics is
Q 2 0 . (1.8)

In TIP the entropy production and the inequality
(1.8) play a far less important role than in the theory
to be presented here. The inequality does serve to place
restrictions on the phenomenological equations which are

obtained according to the following ad hoc procedure:

13

After the energy equation (1.6), the equation of motion
(1.2) and the equations of continuity (1.1) for the fluid
have been introduced into the entropy balance equatiod and
the resulting terms arbitrarily separated into flux and

source terms, the entropy production can be_written as,

r
¢=F:J+Z§-J,, (1.9)

Where terms represented by the symbols F are called forces
and those represented by J are called fluxes. This classi-
fication into fluxes and forces is often arbitrary and it
is impossible to show that any given assignment is correct.

One of the more common assignments is

F = (g + pl) . g = 2}} . (1.10)
F0 = yzn T, J0 = - q, (1.11)
if: u = - -
Ei VT~ , Ji pi(ui u), (1.12)
where
Vu1=Vu-X +§VT (113)
~T i ~ 1 ~i i~ ’ '

and where the heat flux q is related to the total energy

~

flux of equation (1.5) by

r
q = J - DELI + 2J.fi . (1.14)

14

Also, Hi and Si are the partial specific enthalpy and
entrOpy respectively. The linear relations for the fluxes

of (1.10-13) are generally written
(9 + pl) = (¢ - 2/3n) (Y°§)} + 2n sym(Yg) (1.15)

where n is called the shear viscosity and ¢ is called the

bulk or volume viscosity,

r
_ l
- q — (zoown'r + 2 nojvTuj (1.16)
3—1
and
r 1
+ Ji = QiOVEnT + jgl flijvTuj 1=1,...r . (1.17)

The coefficients Qij are called phenomenological or Onsager
coefficients. The number of independent phenomenological
coefficients is limited because of the linear dependence

of the diffusion fluxes
r
2 J. = o (1.18)

and is further restricted by the entropy inequality (1.8).
See Appendix A for a derivation of these restrictions. It
should be made clear that the phenomenological expressions
given here are not unique and that TIP generally considers
many transformations among the forces and fluxes. The pro-

cedure is quite constant however: the phenomenological

15

relations are always determined in terms of those particu-
lar quantities appearing in the expression for the entrOpy
production.

It is certainly obvious that phenomenological
relations must be used in order to make the theory deter-
ministic and that it is convenient to write out the source
terms of the entropy balance equation in order to make use
of the second law inequality (1.8) for the entropy pro-
duction. It is not so obvious, however, that the phenome-
nological relations are to be written only for those quan-
tities called fluxes in the expression for the entropy
production or that they are only determined by linear
combinations of only those forces of the same tensorial
order. For example, TIP suggests that the same expression
used for the stress tensor in a pure fluid is also appro-
priate for a mixture of fluids. Furthermore it is not
obvious that the heat flux should be independent of the
velocities of the components. The motivation for this
ad Egg procedure for obtaining phenomenological eXpressions
in TIP needs to be examined and at the same time the pos-
sibility for making the system deterministic by other
procedures should be considered.

The theory of mixtures to be developed in this_
chapter uses the thermodynamic methods of Coleman and Noll
(1963) and Coleman and Mizel (1963, 1964). This approach

combines the rational use of the entropy inequality (1.8)

16

and the principle of equipresence (An independent variable
present in one phenomenological equation for a material
should be so present in all, unless its presence is in
direct contradiction with the assumed symmetry of the
material or the laws of thermodynamics.) to obtain a sys-
tem of transport equations consistent with ordinary thermo-
statics but without the need for the second part of the
assumption of local equilibrium. This approach to non-
equilibrium thermodynamics, after careful extension to
mixtures, eliminates the inconsistencies and removes most
of the criticisms mentioned above.

Coleman, Noll, and Mizel considered single thermo-
elastic materials, and the extension of their methods has
been attempted by only a few workers. Green and Naghdi
(1965, 1967), Ingram and Erigen (1967), and Bowen (1967)
have made such an extension, but have failed to produce a
theory which is both well motivated physically and also
consistent with classical thermostatics. Stimulated by
the shortcomings of these authors, Mfiller (1968) succeeded
in formulating a thermodynamic theory of mixtures of fluids
which is consistent with both thermostatics and TIP in
those circumstances in which they can be compared. Mfiller's
major contribution was to show that the form of the entropy
production could be deduced from the principles set forth
by Coleman and Noll. The theory of this chapter extends

Mfiller's work in such a way that the comparison to TIP is

17

easily made and so as to make the resolution of the incon-
sistencies mentioned above most obvious.

The development of the theory will be along the
following lines: The conservation laws are introduced in
the standard fashion to obtain the equations of continuity
for each species, the energy balance equation, and in dis-
tinction with TIP, equations of motion for each component
in terms of partial stress tensors. The entropy inequality
and balance equation for entropy are introduced next. The
energy balance, equations of continuity and motion, and
the entropy balance equation are all used to write the
entrOpy inequality in great detail. The (Helmholtz) free
energy is introduced as the energy variable at this stage
of the formulation for convenience in writing the results.
In order to make the information content explicit (with
respect to number of unknowns versus number of equations)
the concept of thermodynamic process is introduced at this
point. It is obvious that the theory at this stage can
not be developed further without phenomenological assump-
tions about the fluid. Having introduced particular con-
stitutive assumptions of an extremely general form (in
particular not linear) the consequences of the entropy
inequality are deduced. One result for this particular
mixture is the establishment of the validity of the thermo-
static functional relations applied in nonequilibrium situ—

ations. Further results are obtained by restricting the

18

class of constitutive relations to those linear in the
independent variables. Such a linear fluid is called an
"ordinary" fluid mixture. After a considerable amount of
algebra and careful application of the entrOpy inequality,
it is found that the equations of transport for the
"ordinary" fluid mixture are surprisingly similar to

those generally found in TIP. There are differences how-
ever; for example, the diffusion fluxes are determined
from equations of motion rather than directly from phe-
nomenological equations. There are some differences due
to the fact that terms involving the kinetic energy of
diffusion are included here but are usually neglected in
TIP (although not necessarily). Finally a complete review
of the problems brought out in this introduction is made
to show that a better understanding of the nature of non-
equilibrium thermodynamics for mixtures has been achieved.

2. Kinematiggand ‘
Equations of Balance1

 

 

In this section the kinematics of motion and the
axioms of balance of mass, linear momentum, moment of
momentum, and energy for a mixture of v components are

considered. Nothing new is presented in this section;

 

lCartesian tensor notation will be employed through-
out this chapter. Vectors are denoted by a single subscript,
vi, while tensors are doubly subscripted, tij. The summa-
tion convention on indices will also be used i.e.,

19

reference could simply be made to the papers of Mfiller
(1968), Bowen (1967), or Truesdell (1957, 1960) were it
not for the fact that notation and some of the well-known
relations need to be introduced to the reader. Accord-
ingly, the derivation of the balance laws is only sketched
here.

Every place in space, xn, is assumed to be occupied
simultaneously by "particles" of all v components. The
density of component a is denoted pa, and v: is the com-
ponent velocity relative to some frame of reference. The

total density of the mixture is defined as2

V
p = {p (2.1)

ov. = X o v9 (2.2)
OI

The diffusion velocity of component u: is defined as

u. = v. - v. (2.3)

 

2For the remainder of this chapter the summation

v
over all species 2 will be denoted 2 .
o=l o

20

The material, substantial, or convected time derivative,

gL , i.e., the time rate of change an observer would detect

if he moved through the mixture with velocity Vi, is re-

lated to the time derivative at a fixed point, §%», by

d_ a a
as- ss” Vi Iii ”-4)

Thegeneral balance equation.--The balance or con-

 

servation laws will be postulated as integral equations.
In regions where the field variables are sufficiently

smooth the integral balance equations may be replaced by
local differential equations. The integral equations of
all the balance laws considered in this chapter are spe-

cial cases of the general balance law

d
EE_]- pawadv = - jr ngsi + Jr paoadV , (2.5)
V S V

where V is a fixed volume with bounding surface S, F: is
the influx of we through S, Go is a volume supply of a
within V, and we can be any component of a tensor, vector,
or scalar. a
8p 31) as. avo.‘
In regions where p w v? o a 0 «——£ ——£ are
o’ a’ 1'5t '5t ’ o'ax.'3x.
continuous, (2.5) may be replaced by the differential

equation

3 8 a
SE 9 W + 3x. (powovj

a _
a C1 3 + fj) _ p O I (206)

21

where the convective contribution to the influx has been

accounted for separately according to

a _ a a

By summing the integral equation of balance for an
individual species over all species the total equation of

balance is obtained,

“‘f
——- pde = - If F.dS. + JI dev , (2.8)

where the following definitions are used:

ow = X puma (2.9)
a
_ a a
Fi _ g (£1 + pawaui) (2.10)
pS = 2 pa a - (2.11)
a
The differential equation associated with (2.8) is
aplp+-§—F=OS . (212)
RE axj j ’ .'

Since appropriate Fj and 8 may always be defined
for any given w, it may be said that all quantities w may
be balanced, by definition. There remains, however, some

indeterminacy in the definitions of the influx and volume

22

source since any arbitrary solenoidal field may be included
in the influx without affecting the balance equation. Even
so the general balance equation is still useful since in
many cases the influx or source may be given a priori, or
special forms for these quantities may be derived from in-
formation about w.

The conservation of mass.--In the absence of chem-

 

ical reactions the mass of an individual component is
assumed to be conserved. Thus the equation of balance

of mass of species a is postulated to be

d _ _ a _

This is a special case of (2.5) with
w=1,f‘?‘=o,o=o. (2.14)

Substitution of (2.14) into (2.6) gives

so
a 3 a _ _
W + K]. {)an — O (a-l,...\)) o (2015)

Summation of (2.15) over all species results in

apv.
Go ;I_
75_t+ 3X _ 0 . (2.16)

23

An alternative form of (2.15) can be obtained by intro-
ducing the rate of change with respect to the barycentric

velocity. Using (2.4), (2.15) may be written in the form

dpa a apa a av?
H'E-I-+ 11). 5357+ p #= 0 (0:1,...V)
J J

. (2.17)

The conservation of linear momentum.--The equation
of balance of linear momentum of species a is postulated

to be
3—.j-pav.:dV'=I-j;(pav:v;.-Sij)dsj -+j;(pm:+pab.:)dV(o=l,...v). (2.18)

This is a special case of (2.5) with

II
I
m
o
o
II

a a
pmi + pobi . (2.19)

aj)= pmi H+p b: (a=l ,...v) (2.20)

The source term here has been split into two parts. The
external body force density pub: gives the part of the
linear momentum source which would be calculated if the
fluid consisted of only component a. The remainder, pmg,
represents a force density due to interaction of component

a with the other components. 89

lj is the partial stress

tensor of component o.

24

Clearly the total linear momentum should not be

affected by internal interactions; thus,
2 m: = o . (2.21)
a

Summation of (2.20) over all species yields by use of

Z pau. = o (2.22)

and (2.21) the following equation expressing the balance
of total linear momentum:
a a

3L
t 3x.
3

j a 1 j i

o. O.

The quantity in brackets corresponds to the total stress

tensor for the fluid as a whole:

a Si. - p u.u. = Z pab. . (2.23)

T.. = Z ( a - pugug) . (2.24)

s..
13 a 13

The conservation of moment of momentum.--The re-

sults needed here are completely embodied in Truesdell's

(1960, p. 546) derivation of Cauchy's second law of motion,

expressed by

When there are no assigned couples and no couple-
stresses, a necessary and sufficient condition for
the balance of moment of momentum in a body where
linear momentum is balanced is

=T.. ; (2.25)

i.e., the stress tensor is symmetric.

25

Hence from (2.24) one concludes that the sum of partial
stresses is also symmetric, although the partial stresses
themselves need not be.

The conservation of total energy.--The balance

equation of total energy is postulated to be

u .ijj

3% f pETdV = - f(q -v..T .+E ij)dS. J+f(pr+2pa VzbaMiV, (2. 25)
V S

where ET is the total energy density, qj is the heat flux
vector, and r is the heat supply density. Equation (2.6)

becomes

_ o
p E- + 5';- (qj - ViTij) - pr 'l' gpavgbi . (2.27)

The internal energy is defined as the difference between
the total and kinetic energy:
E = E - ¥L 2 p v v (2 28)
2p a o o O C

A balance equation for the kinetic energy is obtained by

using equations (2.15), (2.20), (2.22) and (2.16):

d l a a _ 8 o o _ a o
0 5E (53 g pavlvi) — Egj [g vi(S13 pauiuj)
o
l a o o a 1
+ 2 g paululuj] - g 1] ij + p g m u

+ z p Vibe (2.29)
OI

&*4 ,,

26

Subtraction of (2.29) from (2.27) with the use of (2.28)

yields a balance equation for the internal energy:

dEI ahi a a a 3v:
pa-t.—+.5_£3=pr-p§miui+gsijfi—j-, (2.30)

where a new heat flux vector has been defined,

a and l a a a
-pau.uj)+§29 u.u.u. . (2.31)

(S
1 a 1 1 3

330,113

3: The Clausius-
Duhem Inequality_

 

 

The mathematical statement of the second law of
thermodynamics takes the form of the Clausius-Duhem in-
equality or the entropy inequality. The form of the in-
equality used by Coleman and N011 for single materials is
not sufficiently general. Even the various entropy
inequalities prOposed in works of TIP include an influx
of entropy different from that for single materials.
Perhaps some of the shortcomings of the earlier theories
of Green and Naghdi, Ingram and Eringen, and Bowen are'
due to their particular choices for the form of the
entropy production. Mfiller (1967, 1968), in fact, showed
that the form of the entropy production could be derived
explicitly within the domains of the principles set forth

by Coleman and Noll.

27

According to the remarks made in the previous sec-
tion any quantity can be balanced. Thus, the balance of

entropy will be postulated to be

6 f
pndV = - Jf (¢.+pnv.)dS.+.]—

where on is the entropy density, 8 is the temperature, y

gym;— dV , (3.1)

 

 

is the entrepy production due to dissipative effects in V,
and ¢j is the entropy flux through the surface. The main
point in this formulation is that, in contrast to all
previous work, the entrOpy flux is not to be specified
in advance. The entropy flux, sj, is to be regarded as
a quantity for which a constitutive relation has to be
formulated, just as is usually done with the heat flux
and the stress tensor.

From (3.1) and (2.4), the result corresponding to

(2.6) is

dn

—._1 £
9 as" ax.“’°+pY+p ° (3'2)

6
J J

The second law of thermodynamics is made mathe-

matically precise by the following postulate:
y 2 0 . (3.3)

This postulate places restrictions on constitutive equa-

tions of the type to be introduced in the next section.

28

It will be convenient for examining these restrictions to
write the entrOpy production, 7, in detail using the follow-
ing definitions and substitutions:

The (Helmholtz) free energy is defined by

wI = EI - 0n . (3.4)

The difference between the entropy flux for mixtures and
that normally used for single materials (viz., the heat

flux divided by temperature) is defined by

k. = ¢. - —i . (3.5)

Substitution of (3.4) and (3.5) into (3.2) yields, upon

rearrangement,
as ’ 3h. aw 3k. h.
_ I j _ _ I _ d6 j__ j 30
pay - p dt + Bx), pr pat pn dt+eaxj 8 ij ° (3'6)

The equation of balance of internal energy (2.30)

is inserted into (3.6) to give

aw 3k. h. 8v?
3 xj D“ a Xj

The partial stresses can now be decomposed into symmetric

and antisymmetric tensors according to

89. = t8. + :9. (3.8)
l] 13 1]

ta. = i (8“ + 8“.) (3 9)
ij 2 ij j1 '

29

T.. = % (Sq - a )

13 1j ji , . (3.10)

which allows (3.7) to be written

d4) 3k. h
I de 36 a
°9Y=“°ar—‘p”ar+esrl"aisi‘olmju§
3' j a
+ 2 tijd ij + g Tijw ij . (3.11)

In (3.11) the two quantities

3v: av?
(3.12)

_ l
dij _ 2' 5x. + x.
j 1
and

Bug Bug 8p So
a _ l 1 1 B B_ B B

have been introduced with the help of the relation

a .—
2 Ti. — o , (3.14)

which is a result of the symmetry of the total stress
tensor.

Equation (3.11) now gives the entropy production
in a convenient form for use in the entropy inequality

(3.3).

30

4. Thermodynamic Processes

 

By a generalization of the definition of Coleman
and N011 (1963), a thermodynamic process for a mixture is
defined as follows:

A process is a thermodynamic process if it can be

described by a set of (5v + 6) functions

v“ = v°‘ (x t) E - E (x t) ' - (x t)
1 1 1' I ‘ I 1' n " n 1'
ma = ma (x t) 6 - 8 (x t) <1) -<I> (x t)
i. i i’ ' i’ i- j. i’

a - a -

for each a, o=1,...v, which satisfy

a) the balance of mass for each component
(p v.) = 0 (o=l,...v) (2.15)

b) the balance of linear momentum for each component

0.

3%":- 1

a _ o o _
-Sij) — pmi + pobi (o—1,..An (2.20)

c) the balance of angular momentum for the mixture

2 ri‘j = o (3.14)

31

d) the balance of energy for the mixture

dEI ahj a a a 3v:
p E?- + axj = pr - o g mini g Sij 3;; (2.30)

for all xi and t. To specify a thermodynamic process it
suffices to prescribe the (4v + 5) functions pa,v:ml H,S 1j’
,6,hi,n,¢i . The remaining (v + 1) functions bi and r
are then determined from equations (2.20) and (2.30).
In the above definition no reference was made to
the entropy inequality (3.3). Its part in the development
will become clear in section 6.

5. Constitutive Assm mptions
£6} Fluid Mixture

 

 

To illustrate the implications of the theory pro-
posed in sections 2, 3, and 4, it is useful to consider
a special mixture, namely, a mixture of elastic materials
susceptible of diffusion and heat conduction. This case
is sufficiently generaly to illustrate the important fea—
tures of the theory, and it is important since it includes
by further specialization the thermodynamic theory of
single materials and the theory of TIP for mixtures with
diffusion and heat conduction.

The equations of balance and the entropy inequality
suggest that constitutive equations should be prescribed
a o

,t ,T..,m. for each o,(o=l,...v). It is

for wI,n, hi ,k. ij 1] 1

i

assumed that the history of the motions of the components

32

and the history of the temperature within the body con-

sidered determine these quantities as functions of xi and

t. The functional relations that connect these functions

with the histories of motion and temperature are called

constitutive equations; their form characterizes a material.
If, in particular, the constitutive functionals re-

duce to ordinary function of

So
36
pslfiglelfi.lugldgjlw§j (5.1)
J 3
then for the purposes of this chapter, the material will
be called a mixture of fluids.

The constitutive relations are then

- B uY 8 Y
)IJI - ADI (pBIfiI Igi ui (dij twij) (5.2)
= B uY B Y
n n (p8,fi ,0,gi ui 'dij'wij) (5.3)
_ B Y B Y
kj - kj (OB'fi’e’gi’ui’dij’wij) (5.4)
G. _ OI. . f8 uY B (DY =
tij "' tij (DB'f i'e’gi, 'di 1]., (003-) (O. l'ooov) (5.5)
o _ a B B wY _ _
Tij - Tij (p8,fi,6,gi,u ,di lj' %j) (o-l,...v 1) (5.6)
a _ a B Y B wY _ _
mi - mi (pB,fi,e, g J_,ui, ,di 1j’ j) (o—l,...v l) (5.7)
=q (p £889:i quB., 81(3) (58)
qi 1 Bl iI I I I lj' 0
where
Bo
f] = 3x? ’ g] = %E.
J J
and

33

These constitutive equations satisfy the principle
of equipresence, according to which the same independent
variables should appear in all constitutive equations un-
less their presence contradicts the entropy inequality
(3.3) or the assumed symmetry of the material. It should
also be noted that the functions (5.2)-(5.8) depend only
on objective combinations of the velocities and velocity
gradients. This choice was made because the principle of
material objectivity (Noll, 1958) restricts any description
of a material property of a material to one which does not
depend on the frame of reference of the observer, no matter
how he moves.

An admissible thermodynamic process is defined as
a thermodynamic process that is compatible with the con—
stitutive relations (5.2)—(5.8). To every choice of the
component densities pa, component velocities v?, and tem-
perature 6 consistent with the equations of balance of
mass of the components (2.15) and with 6 > 0, there cor-
responds a unique admissible thermodynamic process. For

when pa,vg and 6 are known for all xi and t, then clearly
a
ij

equations then determine wI,n,ki,qi,m:,t:

f:,g3i,u:,d , and ng are determined. The constitutive

a
and Ti Once

J" 3"
these fields are known, then b? and r are determined by
the balance equations.

Thus, at a given time t, it is possible to specify

arbitrarily not only pa,vg, and 6 (subject to (2.15) and

34

de dg dug dd?.
6 > 0) but also the derivatives 3E , dt& , dE_ "—dEl I

dw.. 3g.
1 i .
.331" and gig-at a pOint xi and to be sure that there

exists at least one admissible thermodynamic process cor-

1 It is not necessarily true,

responding to this choice.
however, that for this choice the entropy inequality will

be satisfied. Extending the idea of Coleman and Noll
(1963), it will be required here that the entropy inequality
(3.3) be satisfied for all admissible thermodynamic pro-
cesses. In other words, the constitutive functions are
subjected to the requirement that the entropy production

be non-negative, identically in the independent variables.
The next section will be concerned with deriving necessary
and sufficient conditions that the constitutive functions
must obey in order that the entropy inequality be satis-
fied for all admissible processes.

g. Consequences of the
Entropy Inequality

According to the discussion in the previous sec-
tion, it is necessary to find the restrictions imposed on
the constitutive relations by the entrOpy inequality. The
non-negativity requirement is most easily applied if the
entropy inequality is written in terms of the independent
quantities discussed in the previous section. From the

assumed smoothness of pa,v:, and 6, it follows that

 

1It should be noted that doc and dfg are not in-

dt dt
dependent because of the relations (2.15). This fact will
be used explicitly in the next section.

35

 

 

 

 

 

 

 

 

 

 

 

a
- p ail -Zp (2::d .32 p .igl £1
d dt 8f? dt
1
awI dg. . awI dug
9 (as)% W( T)dti p 3ua '3?-
i
a d.
awI ddii . awI dwl
- Z 0 dt - 2 p ——5 dt (6.1)
a 3d? a 3w..
13' A 13
and
3kg 3kg a 3k2 a 3k
6 Si‘ — 2 e 33— fg + 2 e «at12 + e 56"9z
2 a a 3fi
8k . ak an?
*9 6—2)912+59 ‘3‘ 3x—1
gi a Bug 2
3k2 3d: . 3k 369.
+ 2 e .3 + 2 e ——— ——£l (6.2)
a 3x a 3x
a 3d.. 2 a 3w.. 2
13 13
where for convenience
O. BZDO.
fiz = 6.767.; (6'3)
2
_ a e
912 _ 3x axg (6'4)
' v-1
2 = X (6.5)
a a=l

If now the first and third terms of (3.11) are replaced by
(6.1) and (6.2) respectively, and terms containing the same

independent quantities are collected, then

36

 

 

 

 

 

 

 

 

3w 3w BwI du:
I gI
39 dt and dt1
1
WI 0‘3 ' 34:1 —-1+e3k + er 6
a adOI awII
13 13
where
. 31de -2. 1“: 6(3__..2)_.,]
Sp t a 9 '5— g,
a a a afi
3k 3k 3k an?
+ 2 6 5_&, f: + Z 6 ——£ f“, 'e'——& §_£
a pa 8f? 1 Eu? xfi
i i
a a
3kg adlj . 3kg awlj
+ 6 3x + e a 3x
a aa?. 2 0 3w.. 2
13 13
a a a a a
a a
Since each of the quantities
a a a
d6 dgi dui ddi' dwi. (6 8)

EE ' EE‘ ' dt ' dt ' dt ' 912

can be chosen arbitrarily and independently of any other
term in the inequality, the last with due respect to its
symmetry, the argument given at the end of the last sec-
tion implies that the necessary conditions for non-negative

entrOpy production are:

oeY = oevred 2 0 (6.9)

37

 

 

awI

n = - (55—) (6.10)
BWI

53E._ 0 (6.11)
3¢I

Bui

BWI

— = O (6.13)
(“3)

3w 0 I (6 l4
(XE-j "' o )
3k 8k.

1 i _

L55; 4- 531- - 0 . (6.15)

 

 

 

The inequality (6.9) can be further reduced, but the results
will not be needed until section 7.

Equation (6.10) is a well-known result of thermo-
statics, and its presence here means that it is also valid
in non-equilibrium thermodynamics for a fluid mixture.

Equations (6.1l)-(6.14) imply that the functional

dependence ofdkEmay be written as

(1

WI = LI(pa'e’fi) 0:1,...V (6.16)

In the next section, in which a linear theory is introduced,
the free energy depends only on pa and 6. This together

with some results established by reference to the equilibrium

38

state, yet to be defined, delineates precisely the domains
of validity of the theory of nonequilibrium thermodynamics
derived from the assumption of local equilibrium.

7. Ordinary Binarprluid
Mixtures (The Linear Theory)

Up to this point mixtures of v components have been
considered. However, subsequent formulas would be even
longer and more complicated than they are anyway if more
than two components were considered. Since there is no
conceptual difficulty in extending the theory to a mixture
of any finite number of components, the treatment will be
restricted, without loss of generality, to binary systems.

The development of the theory is further restricted
in this section to the special case of mixtures whose con-

stitutive relations are linear in the variables

Bo 3p
1 2 36 1 1 2 1
FT?" ’5— 73—; ' “2 ' dij ' dij ' wij (7'1)

This special material will be called an ordinary binary
fluid mixture. That such a material is special is obvious
from the restrictions placed on the constitutive relations.
A conviction that this material accurately represents the
binary fluid mixtures commonly encountered by the chemist
in his laboratory is the reason for the label "ordinary."
The appropriateness of this label can only be determined

by testing the ability of the constitutive relations to

39

represent accurately the experimentally measured responses
of the material to variations of the independent variables.
It should be mentioned, however, that under special condi-
tions, according to which it is possible to neglect the
effects of the so called inertial and viscous terms in the
equations of motion, the ordinary binary fluid mixture is
precisely the same material that TIP claims to describe.
This correspondence will be made explicit at the end of
this section.

Material objectivity and linear representation

 

theorems.--So far the principle of material objectivity
has not been used to its full extent; material objectivity
has been used only in the selection of independent vari-
ables. A further consequence of the principle of material
objectivity is that the constitutive functions must be
isotropic scalar, vector, or tensor functions. Repre-
sentation theorems have been worked out for these iso—
tropic functions (see Truesdell and Noll (1965) and Smith
(1965)). These theorems will not be needed here in their
full generality, but for the case of functions linear in
the variables (7.1) the following relations must hold for

scalars s, vectors wi, and tensors tij:

40

B B l _ _
S (DBIfi 16:9 lluill Id lj’ wij) - 5(0816) B - lI2

f8 8 l B l
wi (DB'f 2'9'92'u2'd2j ”(3') = .ngi + 821‘” p6f£+wn “1
2
3 f8 1 B l _ _ B d8
tij (DB! flierggruirdgjrwzj) P51] + le¢ dnnsij
2
+ 2 2n8 {dB }
B=l
a B B l _ _ 1
tij (ps,f£,e,g2,ug,d£j, wzj) — “wij (7.2)

The superscripts s and a on tij indicate the symmetric
and antisymmetric parts of tij respectively. The coef—
ficients WT, wa,wD, p, ¢, n, and 11 may depend on pa and

8; {dij} is the traceless part of d?

ij'

The functional relations given in (7.2) are usually
obtained in TIP by application of Curie's theorem: enti-
ties whose tensorial characters differ by an odd integer
cannot interact in isotropic systems.l Interactions be-
tween quantities of different tensorial orders are well
known in the kinetic theory of gases. The separation

of effects described by the writers on TIP are simply

those that follow by linearization of isotropic functions

 

 

1This theorem has been attributed to P. Curie
(1894), although apparently he neither stated it in this
form nor proved it.

41

and should be considered only within the context of the
linear theory.

Constitutive relations in the linear theory.--In
section 6 some general restrictive conditions on the con-
stitutive relations (5.2)-(5.8) were obtained. These re-
sults are summarized in equations (6.9)-(6.15), and they
will be applied now to the most general possible linear
constitutive equations.

According to the representation theorems (7.2) the
free energy in the linear theory does not depend on density

gradients and thus for an ordinary binary fluid mixture
(I = (I (91' 02. e) (7.3)

Here, use of (6.ll)-(6.l4) has been made in leaving out
the dependence on WI on f?, gi, ug, dgj' and ng. This
means that for an ordinary fluid the free energy has the
same functional dependence as that assumed in TIP on the
basis of the assumption of local equilibrium.

According to (7.2), the heat flux hi in the linear

theory must be a linear function of the vectorial quan-

tities in (7.1). Hence

__ _ 1_ 2_ 1
hi — ATgi lei szi Quj (7.4)

where the coefficients may depend on pl, 02, and 6.

42

Similarly, mi in the linear theory is given by

1 _ _ _ 1 _ 2 _ 1
mi - mogi m fi m fi mDu.

l 2 l , (7.5)

where again the coefficients may depend on p1, p2, and 6.
A linear constitutive relation must also be for-
mulated for the vectorial quantity kj; hence, from (7.2)

and the restriction (6.15) already found for ki, the form

_ l _ l _ 2
ki — + Kui lei lzfi (7.6)

is expected. However, by further rearrangement of the
reduced entropy inequality (this possibility was mentioned
previously) it will be shown later that ki must be inde-
pendent of density gradients also.

Again from (7.2) it follows that the antisymmetric

tensor Tij is proportional to w1.:

lj

l _ _ l
Tij — uwij (7.7)

where u may depend on p1, p2, and 6.
The symmetric tensor tgj is given, according to
(7.2), by

ochB

O.
t..=—p(3..+ leqp dnndij

2
+ 2 2n“8 {dEj} (6:1,2) (7.8)

43

a8 8

Here again the coefficients pa, ¢ , and nu are in general

functions of pl, p2, and 9; {dzj} is the traceless part of

dzj.
Restrictions Imposed on the Linear Constitutive
Relations by the Entropy Inequality.--In order that the
entrOpy inequality be written in terms of independent
quantities, the equations of balance of mass are intro-

duced into the reduced entropy inequality in the form

_ _ a a _ a
aE—‘- uifi padijdij (7.9)

At the same time several of the derivatives which are zero
for the linear case are eliminated and terms involving the
same independent quantities are collected. The resulting

entropy inequality is from (6.7)

 

 

 

 

 

 

 

 

 

 

e = IL 3w: f1 _ 1 361 f2 u1 + 6 3k2)_ fig
9 Yred 9 pl 53; 1 p2 3oz 1 D1 1 56‘ e 91
3k 3k
+ Z 6 g—& f: + 2 6 -—§ f: p(l mi+£ mi plui
a pa a Bfi 2 pl 02
[3k. /3k. 3w
+ e -—% ml + e d). {p 3—5 pad9.6..
Bu. 1] Bu_ 1] a pa 1] 13
l l
_ 2 Bo 8k da - 9 Bkfi £f1-£ f2 0 ul
6 p Bui 13 9 Eu: “5.1 D2 1 l 1
a a l 1 l 1 l
+ Z ijdlj (31 ij + —2le plwij 2 0 (7.10)

44

where the identity

Bug
d l B l B
5;? = “in + dzg— 3 g psdiz - 3 g uEfz (7.11)

has been used.

The quantities fgg are independent of all other
terms in the inequality (taking proper account of their
symmetry) and in the same manner as in section 6 the

following restriction is needed:

Bkz ki
-—a- + —a- = 0 (7.12)
Bfi Bfg

This implies that 21 and 22 are both zero and thus the
most general possible linear constitutive relation for ki
is

k. = Kn: (7.13)

Equilibrium properties.—-Further conclusions may

 

be drawn from the inequality (7.10) without assuming more
particular constitutive equations. These conclusions refer
to the equilibrium state, and they establish a consistency
between this theory and thermostatics.

The material is said to be in equilibrium if the

independent variables

, d.., d.., w.., u?) A=l,...21 (7.14)

45

are all zero. In order to write the entropy inequality
completely in terms of these quantities the linear rela-
tions are inserted into (7.10) to yield

AT Cl 1 C2 2

_ l l
peYred . If'gzgz + TT'ngQ + 7§1g2f£ + p

— + — m p u
01 92) D 1

 

1 1 1
61+ 32)mo]°1“292

 

 

3W
p ( I 6 (3K ) (1 1 ) 6K ] 1 1
+ _ + — —— +;)—1 + — m - —— p u f
[91 501 91 801 91 02 l 091 1 2 2
+ B if; 6 25. + i + l m + 25 ulfz
pol oz 2 002 D1 2 l

 

+ Z 2 ¢“Bd“ d8 + 2 2 2n"‘8 {dgj} {afij}
C1 0.

 

 

8)
1 l u l l I _ l l
+ '61 + ‘2) 61 “’1:- ij +[ o + 901(7)}: P )du
8w er ]
I 2 2 2
+ _ - “— - d 7. 5
[902(3p2) p p 22. ( l )

The entropy production y may thus be considered as
a function of the variables (7.14), and from inspection of

(7.15) it is seen that
Y (X1,X2,...X21) 2 O (7.16)

and

y (0,0,....O) = O . (7.17)

46

Hence y has a minimum value in equilibrium. A necessary

condition for this is

(g?) = o (7.18)
A 0
where the subscript 0 indicates that the function to which
it is attached is to be evaluated in equilibrium.
Application of (7.18)to (7.15) results in the fol-

lowing restrictions:

 

 

Bw GKp
l I 2
P = CD ‘—-'+ (7.19)
l 391 0
8w 6K9
2 I 2
P = 992(§E—)- —Er— (7.20)
If the pressure is defined as
then from (7.19) and (7.20)
aw
I
p = p 2 pa 53— (7.22)
a a

 

 

Eq. (7.22) is the usual definition of pressure in thermo-
statics, and thus provides another example of an equilib—
rium relation that carries over to nonequilibrium processes
in an ordinary binary fluid mixture without need for local

equilibrium assumptions.

47

If the chemical potential is defined as1

 

 

BwI
“a ‘ ”’1 + p “—30 (7.23)
a
then the relation
2 Dana = pr + p (7.24)

O.

is easily obtained, from which the Gibbs-Duhem relation

Zpa—3=-a£.-pn 39- (7.25)

follows by use of (7.3).

The results of the theory given here for the ordinary
binary fluid mixture are thus completely consistent with the
thermostatic theory and all the formalism of the latter can
be taken over for use in nonequilibrium situations.

There are still further consequences of non-
negativeness for entropy production. If, for example,

all variables in (7.14) are chosen to be zero except for

1

wij' then it is necessary that

u 2 O (7.25)

In the same way, if all variables except dgn are chosen

to be zero then

 

1The symbol Ha used here for the chemical potential
should not be confused with the symbol u used previously
for the phenomeological coefficient for the antisymmetric
part of the stress tensor in equation (7.7).

48

12 21
11 22 ll 22 1
¢ 2 0 : ¢ 2 0 , ¢ ¢ - 4 (¢ +¢ )3 0 (7.26)
One also finds using the same technique
U11 2 0 ' r12 2 0 ' r111r122 __ %_ (n12+n21)3 o (7.27)

If now the temperature gradient 92 is chosen to be some

arbitrary vector Xi' f: and f:

to some other arbitrary vector Yi’ and all other independent

are both chosen to be equal

variables are taken to be zero, then it can be seen by com-
pleting the square that in order for the entropy production
to remain non-negative for all choices of Xi and Yi it is

necessary that
A 2 0 (7.28)

and

C + C = 0 (7.29)

Furthermore, if the density gradient ft is chosen to be
equal to some arbitrary vector Xi' gfl and f: are both
chosen to be equal to some other arbitrary vector Yi' and
all other independent variables are taken to be zero, then

by completing the square, it is necessary that
C = 0 . (7.30)
Eqs. (7.29) and (7.30) imply further

C = O . (7.31)

49

Precisely the same technique used with the diffusive velocity
ufi, instead of the temperature gradient, and the two density

gradients yields the three necessary conditions

 

2

p 8 8K 1 1 8K
——— m + — + ——-p - —— = o (7.32)
0192 1 91 591 02 Q2

1 1

_E_z_m+§.§§.-.}_p2=o (733)
p192 2 p1 3p2 p3
mD 2 o (7.34)

Finally if all the independent variables except for the

temperature gradient 92 and diffusive velocity u:

taken to be zero, the non-negativeness of the entropy pro-

are
duction requires

9 meD 1 8K Q p z 0 . (7.35)
—T‘z 9— 6"

Summary of the linear constitutive relations.--A

 

summary of the most general linear relations possible for

an ordinary binary fluid mixture are

h = - x g - Qul (7 36)
i T i i '
l _ _ 1 _ 2 _ 1

mi — mOgi mlfi m2fi mDui (7.5)

where

11

ll
0

\V

\V

\V

\V

\V

Ku;
1

 

 

 

 

U ij
- paé + E oaBdB 6.
ij B=
a=l,2
801 GKp2
+
001 g3; p
BwI er2
pp '55;- o
9.11-2 (9.15.)-
pi 01 3pl
—.e_(3K+._l_.p
91532 03
0
0
1 3K Q
71'[6T)+6+
0
22

\V

(7.13)

(7.7)

(7.8)

(7.19)

(7.20)

(7.32)

(7.33)

(7.28)

(7.34)

(7.35)

(7.25)

(7.26)

(7.27)

51

These relations are necessary and sufficient for the entrOpy
production always to remain non-negative for any admissible

process in an ordinary binary fluid mixture.

8. The Transport Equations

This section is intended as a summary of the dif-
ferential equations governing the evolution of the state
of the ordinary binary fluid mixture in thermodynamic
processes. No new results are derived here. The equa-
tions are rearranged and expressed in a form more convenient
for applications. This includes replacing p1 and p2 in
favor of W1 and O as independent variables, where the mass

fraction w1 is given by

wl = 01/0 (8.1)

The balance of mass.--The two independent equations

 

of conservation of mass to be used here are

do _ _ l
HIE" p — (8.2)

which follows from (2.16) upon use of (2.4), and
.1
dw Bji

__1==_ ___
p t 8x.

d (8.3)
1

which follows from (2.15) upon use of (2.4) and the defini—

tion of the diffusion flux

jj== u) (8.4)

52

The balance of momentum.--The overall equation of

 

motion of the mixture is derived from equation (2.23):

 

dv p(V'-V ) 3w 3n..
1 _ a _ B 36 _ 1 8p _ l 2 l ij
p at - 2 pubi. 8' 5x. 98' 5x. 8' 8x. + x. (8‘5)
a 1 1 1 j

and the following thermodynamic identity has been used for
the pressure gradient:

_§£ ____ B 88_1 a _ p(Vl-VZ) awl (8 7)
3x. FT'fiil 68‘ x. 8' 3x. ' °

1. l l J.

 

where 8 is the thermal expansivity,

B=—.];§.B.
9 3T
le1
8' is the isothermal compressibility,
811%.
p 3p T,w

 

and Va is the partial specific volume.
The equation of balance of linear momentum (2.20)

is written more conveniently

(8.6)

(8.8)

(8.9)

613'1
p i _ D .1 38 3
- —— 1+ D §——+— (u u)
l 8 l l l 3 2 2
-[§;3—x;(s +P513"6;5§;‘51j+1° 13’]

 

 

J 3
9192 3"j p2 5"j
)
+ .1.1 1 8p .1 (9(92'91)+°192 awi); (8 10)
Jijj EIEE'EET’ 9\ 2 2 8§.y '
where the coefficients D1 and D2 are related to mo, mD,
and K by
psz
D1 = 2 (8.11)
D192
2
p m
_ 0 K 8 8K — _-—
D2 - T— + 6- + B—(fi <8]. 52) (8.12)
0102 1 l

The partial specific entropies, 85, have been used in (8.12).
The special derivative of chemical potential used in (8.10)

is defined by

 

 

— _ .— -— 2
E3(V1"’2) 88 +[0(V1‘V2’ pull] 8“1
l

 

3 _
axi (“l-“2)T _ 8' 3x. 8' + p2 Exi
+ (Vl-VZ) 39 (8 13)
08' 3x. °
1
where
U = iii (8 14)
11 8x1 8,p

 

 

54

The heat flux.--A new heat flux is introduced by

 

the definition

— (fi'l - HZ) j? (8.15)

Q- = Q: J

3

where Ha is partial specific enthalpy defined by

- _ — l a
H - 65a + “a + f uiu. (8.16)

a
0'. 1

From the definition (2.31) of the heat flux hi and the

constitutive relation (7.36), qj is seen to be

1 _ 80 .l a a a-
- qJ - Bl axj + BZJj + g ui (Sij + p Gij) (8.17)
where the coefficients B2 are related to AT and Q by
B1 = IT (8.18)
and
B = 0/8 + 6(5' - §') - 95 (8 19)
2 1 l 2 p:L I '

The temperature equation.--In section 5 it was

 

determined that a thermodynamic process could be specified

by the fields wl, p, ji, Vi' and 6. So far differential

equations have been given for all of these except 8. The

temperature equation can be obtained from the following

thermodynamic identity:

— de_ 211 8.8-9.3- [—_— _
QCV 376“ 99 dt+p8' dt (38(81 52) 8'

 

55

The equations of balance of mass (8.2) and (8.3) will be
substituted into the last two terms of (8.23) and the first

term will be replaced by the entropy balance

dn _ _ 8
p6 3E — 5E3'6¢j + @jgj + pay + or (3.2)
The entrOpy flux to be used in (3.2) is determined
from the definition (3.5) of ki together with the constitu-
tive relation for ki, the definition (2.31) of hi and the

definition (8.18) of qi:

_ l — _— .1 a a a
e¢j — qj + 8(sl $2) jj + g ui (sij + p Gij) (8.21)
Similarly, from equation (7.10) for the entropy production,
the relations (2.20), (7.12), (7.23), (7.19), and (7.20)
may be used to write the entropy production in the more

convenient form:

 

 

_ 8 .1 a a a a
NY 7 ‘ 45-93- [5:]— (“1 112)] 3j " g ”‘1 65; “1
§(Sij + p dij) 5;; (8.22)

Substitution of these last two results into (3.2) gives

dt 8x 3 1 2) 8x. 33 + or
J 3
8v.
.1 3 —v —v 1
—Jj §;—-(Hl - H2) + Fl] 3X (8°23)

56

where the system has been restricted to external fields

derivable from a potential

a 39a
bi = - 5?- (8.24)
1
and where
—I _
a - Ha + “a (8.25)

The final temperature equation is now completed by
inserting (8.26), (8.2) and (8.3) into equation (8.23).

The result is

 

5 d8 = _ a _ 65(V1‘V2) a .1 + r
p v a’E _3xj qj 8' ij Jj 9
av. av.
86 1 _ .1 3 ' _ ' 1
' 87'52; 31 6§;"H1 Hz’ + "ij IE; (8'26)

Summary of transport equations.--The complete set

 

of equations which are necessary to describe any admissible
thermodynamic process in an ordinary binary fluid mixture
are equations (8.2), (8.3), (8.5), (8.6), (8.10), (8.17),

and (8.26) together with the apprOpriate boundary conditions.

9. Discussion

 

It is not surprising that the transport equations
for the ordinary fluid mixture bear striking resemblances
to the transport equations of TIP (see Horne, 1966) since
hopefully they describe the same material. Briefly, a

comparison of the equations yields: (1) The equations of

57

continuity of mass, (8.2-3), are identical with those of
TIP; (2) The temperature equation (8.23) (derived from

the energy balance equation) is identical in form to that
of TIP except that the possibility of external sources of
radiation has been included; (3) The phenomenological
equation for the heat flux (8.17).contains viscous terms
and a kinetic energy of diffusion not found in TIP; (4)

The overall equation of motion (8.5) is identical in form
to that of TIP; (5) The viscous pressure tensor (8.6)
appearing in the equation of motion contains a contribu—
tion from the kinetic energy of diffusion; (7) The phenom-
enology of the partial stress tensors (7.7-8) provides
completely new expressions for the determination of stress;
and (8) The equation for the diffusion flux (8.10), usually
given as a phenomenological equation in TIP, is determined
here from the componentwise equations of motion. If the
inertial and viscous terms of (8.10) are neglected, the
equation is identical with that of TIP.

On the surface it appears that the major differ-
ences between the two formulations arise because of the
use of the partial stress tensors and componentwise equa-
tions of motion. It should be apparent, however, that
many of the underlying inconsistencies have also been
eliminated. For example, the careful application of Cole—
man's thermodynamic methods has eliminated the need for

the local equilibrium assumption as it is commonly used.

58

An extension of Mfiller's work has made it possible to
determine the entrOpy flux by derivation rather than by
making an arbitrary separation of flux and source terms

as is done in TIP. It has been shown that the entropy
production may be written as a quadratic form which allows
the second law entropy inequality to be used with maximum
effectiveness. Finally, the equivalence of the pressure
appearing in the Gibbisian equation for the energy and

the negative trace of the stress tensor has been estab-
lished.

One result of the identification of pressure and
the ability to write the entropy production in a quadratic
form is that is is now possible, for the first time to
give an unequivocal criterion for mechanical equilibrium:

A system is in mechanical equilibrium if and only if

These restrictions also imply that at mechanical equilibrium
the external forces are completely balanced by the pressure

gradient in the fluid.

CHAPTER III

PERTURBATION EXPANSION METHODS

1. Introduction

 

The particular stimulus for the research in this
thesis has been the need for a better understanding and
interpretation of the results of pure thermal diffusion
experiments. The previous chapter is the culmination of
an effort to better understand the fundamental theory of
irreversible processes in liquid mixtures. The following
quote from Bellman (1964, p. 1) sets the tone for this
chapter:

. . ., one of the major problems confronting the
mathematician, after he has passed the first hurdle of
achieving a reasonable analytic formulation of a physi-
cal process, is that of deriving useful and meaningful
approximations to the solutions of the equations de-
scribing the process. In some cases, analytic ingenu—
ity, alone or abetted by digital computers, will furnish
the desired eXpressions; in other cases, a combination
of "low cunning" and physical intuition will provide
the essential key.

In still other cases, a completely new interpreta-
tion and formulation of the physical situation is re-
quired. A major obstacle at the start of any research
is ignorance of where the real difficulties lie. The
greater part of the time, only perseverance and plodding
effort yield this vital information.

Indeed, there is no reason to believe that the research

reported here has been any different.

59

60

In this section, prior to a discussion of particular
problems, the simplest and most useful of all approximation
techniques will be presented: the expansion of a solution
as a power series in a parameter--the classical perturba—
tion technique upon which much of the edifice of science
rests.

The basic idea of the perturbation technique can
be exhibited most easily in abstract terms. The introduc—
tion by Bellman (1964) is probably the best starting point:

Suppose that we are required to solve the equation

N(u) = v (1.1)

which for any of a variety of reasons, is inconvenient to
tackle in its original form. It may, for example be non-
linear, be linear but of high dimension, or possess variable
coefficients.

Let

L(u) = V (1.2)

be an auxiliary equation that possesses a useful explicit

solution

u = T(v) (1.3)

in general, the unique solution of (1.2). In practice,
this means that L(u) is a linear operation on u. Let us

henceforth assume that L is linear.

61

Returning to (1.1), we write it in the form

L(u) = V + R(u) (1.4)
where the function R has been defined as

R(u) = L(u) - N(u). (1.5)

To facilitate our discussion of particular solutions
of (1.4), we introduce a parameter s and consider the new

equation
L(u) = v + eR(u) (1.6)

In some situations, the introduction of e is solely
a mathematical artifice that permits us to do various
types of "bookkeeping" in a systematic fashion. For
example, it allows us to group terms of comparable degrees
of approximation in a methodical and convenient fashion.
In a large number of situations, however, this parameter
occurs naturally, representing some physical quantity.
In the application to pure thermal diffusion experiments
we shall have use for both interpretations.

It is quite natural and sensible to begin with the
study of those equations where s is close to zero. Let us

then look for a solution of (1.6) having the form

u = u + eu + ezu + ..., (1.7)

a power series in s with coefficients that are independent

62

of 8. Clearly, the leading term uo, obtained by setting

6 = 0, is a solution of the linear approximation
L(uo) = v (1.8)

a solution that we shall write as u0 = T(v).
To obtain the subsequent coefficients, ul, uz, ...,
in a systematic fashion, we substitute u as given-by (1.7)

into (1.6) and equate terms, obtaining thereby

+ ezu + ...) = v + eR(u0 + an + ...) (1.9)

L(u0 + eul 2 1

Since L is by assumption a linear Operator, the left side

becomes
2
L(uo) + sL(u1) + e L(uz) + ... (1.10)

Assuming that R(u) is analytic in u so that we can expand

the right side of (1.9) as a power series in e, we have

2 ._
R(u0 + eul + e u2 + ...) — R(uo) + eRl(u0,ul)

2
+ e R2(u0,ul,u2)+... (1.11)

where, as indicated, the coefficient of uk depends only
upon the quantities un, n s k. Combining the expansions of
(1.10) and (1.11), and equating coefficients of e, the
single equation of (1.6) gives rise to the infinite system

of equations

63

L(uo) = v
L(ul) = R(uo)
L(uz) = Rl(uo, ul) (1.12)

L(“1;“) = Rk(“o' “1' °'°' uk)

and so on.

The important point to observe is that this system
of equations can be solved recursively; that is, the deter-
mination of uk involves a knowledge of un, n c k.

The first equation yields

u0 = T(v). (1.13)

From the second, we derive the relation

u = T(R(u0)) = T(R(T(v))). (1.14)

1

Continuing in this way, we see that we can express each

u solely in terms of v.

k
The infinite series in (1.7), whose coefficients are

determined in the foregoing fashion is called a formal solu-

tion of the original equation in (1.6). To obtain a formal

solution of (1.4), we need only set 6 = 1.

Pure thermal diffusion experiments.--Historically,

the application of nonequilibrium thermodynamics to the

problem of interpreting pure thermal diffusion experiments

64

has been made by considering the magnitude and variability
of the terms in the transport equations. Among the assump-
tions that have been made in formulating working equations
in the past are: (l) constant phenomenological coefficients;
(2) a one-dimensional temperature gradient; (3) a density
independent of composition; (4) no convective motion; (5)
no viscous effects; and (6) dilute solutions. Undoubtedly
there is some validity to these approximations in some
cases, and their use sometimes provides an adequate first
order theory. However, it has been impossible to provide
corrections to the theory for the cases in which there is
significant deviation from these assumptions and it has
thus been impossible to interpret properly the experimental
results. The use of the perturbation technique not only
provides these corrections but at the same time facilitates
the solution to the transport equations.

It will be seen from the description of the pure
thermal diffusion experiment in the next section that two
characteristics stand out as natural possibilities for the
application of the perturbation method: (1) The phenome-
nological coefficients are fairly constant and any varia-
tion might be taken as a perturbation. (2) No steady state
convection is expected if it is possible to minimize heat
losses; convective corrections are therefore introduced by
the use of a heat loss parameter. In fact, as will be seen
in section 3, we employ a double perturbation expansion to

account for both effects simultaneously.

65

2. The Experiment

 

The experimental device used in studies of pure
thermal diffusion is usually described as a sandwich-type
cell consisting of two thermostated, horizontal metal
plates with an arrangement for containing a homogeneous
binary liquid solution of known initial composition between
the plates. A temperature difference is maintained between
the plates by the external thermostatting mechanism with
the warmer plate situated above in order to minimize con-
vective effects. The temperature gradients which are thus
present in the liquid result in more than just heat con-
duction. The two components in the mixture will tend to
separate or demix under the influence of the temperature
gradient (thermal diffusion). Of course the resulting con-
centration gradient produced by this effect immediately
induces ordinary diffusion which tends to remix the com-
ponents. The total diffusion flux will be observed as the
resultant sum of these two Opposing fluxes. The funda-
mental measurement in the pure thermal diffusion experiment
is the composition distribution in the cell under steady
state conditions (all local time derivatives are zero) for
which the thermal diffusion and ordinary diffusion fluxes
are as closely balanced as possible (total diffusion flux
approaching zero) and convective effects are at a minimum.

The boundary conditions.--Ideally the boundary con-

 

ditions to be used with the temperature equation would

66

consist of the temperatures at the top and bottom of the
cell, together with a condition that the heat flux at the
side walls is zero. Practically these boundary conditions
cannot always be obtained and any observations of convec-
tion in the cell are probably for this reason. For the
purpose of this treatment it will be assumed that the
temperatures of the top and bottom plate can be controlled
as accurately as desired, but that there is a slight heat
loss at the side walls. There are two principal reasons
for the constant temperature choice. First, the primary
objective of this chapter is to show how a perturbation
theory can be formulated for the heat loss problem, and
the variable temperature problem only complicates the
situation. Second, there are many practical experiments
that can be done for which the constant temperature assump—
tion is quite good. In any case if the inclusion of vari—
able temperature at the horizontal surfaces is required,
the procedure will not be very much different from that
which will be introduced for handling the heat loss prob-
lem. With these considerations in mind the boundary con-

ditions for the temperature equation are:

8 = Tm + 1/2 (AT) at z = h (2.1)

6 = Tm - 1/2 (AT) at z = 0 (2.2)
and

§_8_=€(T —8) atr=a (2.3)

67

where it has been assumed that the z axis is vertical, with
the height of the cell given by h, and that the fluid is
contained between the plates by a cylindrical surface of
radius a. The difference between the temperature at the
top and at the bottom has been denoted by (AT) while the
average of these two temperatures has been denoted by Tm'
The heat loss coefficient is given above as e, and the sys-
tem is assumed to be radiating heat into surroundings which
are at a temperature of Tr’ Finally, we take it as an
assumption that all of the prOperties of the fluid are
axially symmetric; 1121' they depend only on the vertical
position in the cell and on the distance from the center.

Although s is a parameter which depends on the
apparatus and must bedetermined experimentally in each
individual case, it is helpful to have an order of magni-
tude estimate. If the fluid is assumed to be contained
within glass walls (thermal conductivity = 0.0028 cal/(sec)
(cm2)(°C/cm)) and to lose heat through a stagnant layer of
air (thermal conductivity = 0.000053 cal/(sec)(cm2)(°C/cm))
of thickness 0.5 cm, then s is given by (0.000053)/(0.5 -
0.0028) = 0.038cm‘l. (See Carslaw and Jaeger (1959) for
more details concerning the use of heat loss boundary con—
ditions.)

In the absence of chemical reactions, the boundary

conditions for density and mass fraction are simply state-

ments of the conservation of mass:

68

Zn h a
j, j” .J- p(r,z) rdrdzdw = M (2.4)
0 0 0
and
2“ h a
‘[ —[ -f- p(r,z) W1 (r,z) rdrdzdw = M1 (2.5)
0 0 0

where M is the total mass contained in the cell and M1 is
the total mass of component one. These two masses are
easily obtained by weighing or density determinations at
the beginning of the experiment.

Boundary conditions for the velocity and diffusion
flux are easily written to insure that no material either
flows or diffuses through a rigid wall. These conditions

are

2 jlz = 0 at z = o, h (2.6)

ur = jlr = 0 at r = a (2.7)

It will be further assumed that the fluid sticks to all
rigid walls. This is the usual assumption in all hydro-
dynamic treatments. This leads to the further conditions
ur = 0 at z = 0, h (2.8)
u2 = 0 at r = a (2.9)
Under certain circumstances it might be possible
that the number of boundary conditions is insufficient in
order to determine properly all the solutions to the dif—

ferential equations. If these circumstances arise, it will

always be possible to ask the experimenter to provide the

69

additional measurements of concentration or concentration
gradient at one or more points in the cell.

Pure Thermal Diffusion

3.
Tansport quatibnsl

 

The equations to be used for the description of the
pure thermal diffusion experiment are those listed at the
end of Chapter II. There are a few modifications and
specializations of these equations which make them more
appropriate for pure thermal diffusion. The equations are
written here in such a way as to make the comparison with
the standard treatments of nonequilibrium thermodynamics
convenient. In particular this involves the following
major assumptions: (1) The viscous terms in both the
phenomenological equation for the heat flux and the equa-
tion of motion which gives the diffusion flux are neglected.
(2) The viscous stress tensor is assumed to be given by
the Newtonian formula; 1121! the stress depends only on
the convective motion. The following simplifications have
also been included: (1) Since the description is to be of
asfieady phenomena, all local time derivatives are set equal
to zero. (2) The only external force to be considered is
gravity. (3) It is assumed that there are no external

sources of radiation.

 

1For the remainder of this chapter direct, as dis-
tinguished from Cartesian or component, tensor notation
will be used. Vectors are denoted v, tensors by t, and

the gradient operator by V = 8

ml

70

The steady state equations of transport, modified

as above, are

v ° Vp + pV - v = 0

~

 

 

_JL_ . - ° _
p p Y Y21 + D1.3.1 + YT(“1
1 2
[p(vl-V2)gh pullo]
+ _
AT p2

 

pEvv ° VG + V . q' + 23 V - v - 2 :VV
88(V'-7')
l 2 .
+ BT Y 0 21 +
_ ' = °
3 3139 + B2V21

(1:)

Notice that the phenomenological coefficient D

pv ° Vv + Vp - pgA - V -'n= 0

“2)

V6 = 0

= (8 - 2/3n)(y '9) l + Zn sym (Y9)

(3.1)

(3.2)

(3.3)

(3.4)

(3.5)

(3.6)

(3.7)

of equation (II.8.10) has been eliminated in favor of the

coefficient 0 which is defined by

_ "2(Vi’vz)gh 02132
(AT)

 

“11 pu11

(3.8)

71

The coefficient 0 will be called the thermal diffusion
parameter. It is related to the thermal diffusion factor

a used in many nonequilibrium thermodynamic studies of

1
thermal diffusion. For example, from the phenomenological

equations given by Horne and Bearman (1967) it is easy to

 

see that
wlwzal = O _ "2(V1"V2’gh (3 9)
5 011(AT) °

The last term in this equation is typically very much
smaller than 0, and the behaviors of 0 and al are there-
fore very strongly correlated. Notice also that in equation
(3.7) the_common coefficient of shear viscosity, n, and

the coefficient of bulk or volume viscosity, 0, have been
used. Finally, in the equation of motion (3.4) the exter-
nal gravitational field has been assumed to be of magnitude
g and to act in the direction 1. In all that follows, it

is assumed that this direction is vertical and
1.: (0,0, - 1) (3.10)

For completeness, the thermodynamic expansions for
the pressure gradient and chemical potential gradient
needed in equations (3.3) and (3.4) are repeated from the

previous chapter:

Vp=-B-TY8+—-VQ+-——ET——le (3.11)

Vw (3.12)

 

4. The Perturbation Expansions

 

In this section the beginnings of the perturbation
method for the pure thermal diffusion method are developed.
In the introduction to this chapter it was suggested that
in a typical pure thermal diffusion experiment the diffusion
flux is small, that if there is only a small heat loss the
convective velocity is small, and finally that the vari-
ability of the phenomenological coefficients is very slight.
While it might be true that the phenomenological coeffi-
cients are fairly constant in experiments in which the
gradients of temperature, pressure and composition are
small, the neglect of their variation in the transport
equations constitutes a serious introduction of systematic
error in many cases. Of course the differential equations
are much easier to solve if they have constant coefficients
and not infrequently it might be suspected that the cri-
terion of constancy is more strongly influenced by the re—
duction of work in solving the equations rather than by
any physical evidence of constancy.

In any case, the perturbation method is particularly
well suited to handle problems of this type. Its use

focuses attention on the major processes while at the same

73

time making possible corrections for the minor effects of
variable coefficients. This is accomplished by expanding
each of the coefficients in a thermodynamic series in the
manner of Hurle, Mullin, and Pike (1965) and introducing
the parameter y for the purpose of ordering the terms.

The typical coefficient is expanded in the following man-

ner:
L = L0 + y[LT(6-60)+Lc(wl-wlo)+Lp(0'00)] + 0(y2) (4.1)
Where
L0 = L(GO'WIO'DO) (4.2)
LT = 8%) 0=00,w1=w10,p=p0 (4'3)
LC (381') e0'w10'90 (4.4)
and
L0 = (3%) eo'wio'po (4.5)

Corresponding to the ordering of terms by the parameter y,
each of the independent variables is given by a power series
expansion in y, so that the contribution of each term in
the series to the total is proportional to the effects of
the perturbation at that order.

It is easy to see that if the differential equa-

tions are solved with the parameter Y set equal to zero,

74

then the solution is that for constant coefficients, whereas
if Y equals one, then the full variability of the coeffi-
cients is taken into account. This same analysis can be
applied to the heat loss a as a parameter. If all the
independent variables are eXpanded in a power series in s,
then the solutions to the differential equations with 6
equal to zero correspond to the zero convection problem.

The perturbation expansions for the independent variables
thus become double power series expansions in Y and e of

the following form:

00 oo.k
8 = T + 2 1 le T. (4.6)
r i=0 k=0 1k
M °° °°ik
p = —7—— + 2 2 6 Y d. (4.7)
a uh i=0 k=0 1k
M 00 oo
1 i k
w = ——- + 2 2 e Y C. (4.8)
1 M i=0 k=0 1k
. co mik
j = Z Z 8 Y J. (4.9)
”1 i=1 k=0 ”1k
” m i k
v = Z X 8 Y v. (4.10)
” i=1 k=0 ”1k

where the first term on the right in (4.6), (4.7), and
(4.8) has been introduced for the purpose of making the
boundary conditions homogeneous. The summations on i in
(4.9) and (4.10) begin with one in order to indicate that
in the case of zero heat loss no convective or diffusive

effects are expected.

expanded

5.

begins with the transport equations listed in section 3.

The expansions (4.6)

75

according to (4.1) in the following way:

0' +

80 +

86 +

Y[OT(6-TO)+0C(wl-wlo)+op(0-00)]
YIDT(6-T0)+Dc(wl-wlo)+Dp(0-00)]
YlKT(6-T0)+Kc(wl-w10)+Kp(0-00)]
YIQT(9-TO)+QC(w1-wlo)+op(0-00)]
Y[VT(6’T0)+Vc(w1-w10)+vp(9-90)1
Y[HT(6-To)+HC(wl-wlo)+Hp(o-po)]

Y[BT(6-RO)+BC(wl-wlo)+8p(0-00)]

+

+

+

+

+

The phenomenological coefficients to be used are

0(Y2)
0(Y2)
0(Y2)
0(Y2)
0(Y2)
0(Y2)

0(72)

Y[8&(e-To)+sé(wl-wlo)+85(p-po)1 + 0(y2)

Cv0 + Y[CVT(6‘T0)+Cvc(w1‘w10’+CVp(9‘90)1

+ 0(Y2)

“0 + Y[0T(8-TO)+uC(wl-wlo)+up(0-00)] + 0(Y2)

no + Y[0T(6-T0)+nc(w1-wlo)+np(0-00)] + 0(y2)

(0 + Y[¢T‘9'To’+¢c‘wl‘wlo)+¢p‘p‘po’1 + 0(y2)

The Perturbation Equations

(4.11)

(4.12)

(4.13)

(4.14)

(4.15)

(4.16)

(4.17)

(4.18)

(4.19)

(4.20)

(4.21)

(4.22)

The analysis of the pure thermal diffusion experiment

to (4.22) are substituted into each of

the equations in section 3 and into the boundary conditions

of section 2.

Terms of the same order in6:and~onr each i and

76

k are collected and separated. The equations at each order
may then be solved independently of the solutions of any
higher order.

The zeroth-order equations.--After the substitu-
tions indicated above, the terms of the equations contain-

ing neither parameter are rearranged to give

 

 

 

2 _
v T00 — 0 (5.1)
th
- _ o A_T,
Ycoo " UOYToo ‘(ZTT “—0. [YToo + (11h) (5'2)
_ M I _
Ydoo — d00 + —§——)[8091 (V000+80)YT00]
a Nb
2
V gh
0 M AT
+ "TH T doo + "2— [YToo + (Ta—)5] (5'3)
0 a nh
The boundary conditions become
T = T — T + l (AT) at z = h (5 4)
00 m r 2 °
T = T - T - i (AT) at z = 0 (5 5)
00 m r 2 '
3T
00 _ _
3r — 0 at r — a (5.6)

for temperature, and

Zn h a

-I- .f .[ doordrdzdw = 0 (5.7)

O O 0

77

C

 

rdrdzdw = 0 (5.8)

 

for composition and density.

The zeroth-order solutions.--The zeroth-order

 

equations have been arranged so that the composition and
density distributions may be easily found from the tempera-
ture distribution. The equation for temperature to zeroth-
order (5.1) has a general solution consisting of linear
terms together with periodic terms appropriate to the
boundary conditions (5.4-6). The nature of these boundary
conditions eliminates the possibility of periodic solutions
to (5.1). The temperature can therefore be given by a

linear function which according to (5.4-6) is

TOO = 1%21 (z - h/2) + (Tm - Tr) (5.9)

The gradient of the zeroth-order temperature solution is

most conveniently written

YToo = ' (%§)§ (5'10)

The density distribution is next determined from
equation (5.3). The solution is easily seen to be an
exponential which satisfies the boundary condition (5.7).

This solution is

78

 

d = M M exP [A(z - h/2)] - 1 (5.11)
00 .72.. 2 .....438.
where
_ a _ A:
A — 80g (Vooo+ BO)(h.) (5.12)

Finally, if the temperature solution is substituted
into equation (5.2) and the density solution used in the
boundary condition (5.8) the following solution can be

found for the composition distribution:

(2 - h/2) + % - — hcoth(2) (5.13)

00(AT)
C =—
00 h
In most practical problems the quantity A has a

-2

magnitude of about 10 , and the hyperbolic functions can

conveniently be replaced by their expansions:

sinh(%P) = (fig) + %(§E)3 ... (5.14)
coth (1:31) (‘33—)-l + 1(5h—)- 1 (53f (5.15)

In using these expressions the cubic terms (as well as all
higher terms) are negligible with respect to the first.
Insertion of (5.14) into (5.12) and (5.15) into (5.13)

yields

 

800 = 2M exp [A(z - h/2)] - 1 (5.16)
a fih

and
00(AT)

2
_ Ah .

79

after the higher order terms of the eXpansions have been
dropped. It must be remembered that these are only approx-
imate expressions for use in practical calculations where
the magnitude of A is small.

Equations first order in Y.--After the perturba-
tion expansions have been introduced into the transport
equations the terms of these equations which multiply Y
represent the corrections to the basic processes due to
variability of phenomenological coefficients. The intro-
duction of the perturbation expansions (4.6) to (4.10)
into the thermodynamic expansions (4.1) of the phenomeno—

logical coefficients gives for the typical coefficient
2
L = LO + YL1 + 0(Y ) (5.18)

where Ll has been introduced for compactness of notation

and is defined as

 

M
_ _ _l...
L1 ‘ LT(TOO + Tr To) + LC (coo + M w10
+ L d + M — p (5 19)
p 00 2 0 °
a Uh

The transport equations of first order in Y are

2 -(A—T) cm (5.20)

1__
01 KO

__1
dz
VCOl == 0 lVTO m[— AT 930 YTOl (5.21)

80

 

 

2
M 8' v0gh
Yd01 ‘ ‘ Adelé + doo + 22;; 195 ‘ V000 + Bo ' EETKTT'YT01
- (v100 + voo1 + 81) YTOO (5.22)

From (2.1) to (2.8) the corresponding boundary conditions

are

T01 = 0 at z = O (5.23)
Tol = 0 at z = h (5.24)
3'1.01 = 0 at r = a (5.25)
3r

for temperature and

 

2w h a

J( ,[ df~ delrdrdZdw = 0 (5.26)
0 0 0

2H h a

J[ _[- _[. Cood01 + doo + 34 C01 rdrdzdw=-o (5.27)
0 0 o a nh

 

 

for density and composition.

Solutions first order in y.--In the same way as in

 

the zeroth-order case the boundary conditions on tempera-
ture admit no periodic solutions to (5.20). Under these

conditions the only possible solution for the temperature

 

is
K 2
Tol =(%§)K§ ggig 2 -(%§) 3%; {KT + KCOO} z(z—h)
_A%
+(93) KOMe (1 - eAZ) (5 28)

2 . h
2K0a n51nh A?

81

while the gradient is most conveniently written

K 2
YT01= " (éf1$)fi 221E; -‘él712 1T1; {KT + KCOO} (2 ’ g)
h
-A—
%) K pMAe Zih eAz l (5.29)

 

2Ko a 2nsinh A7

Substitution of this expression for the temperature
gradient into equation (5.22) for the density together with
the detailed expressions for the phenomenological coeffi-

cients, allows equation (5.22) to be written

 

dd
01 _ AZ M
dz ... Adol + [B + 2C2 + ADe ](d00 + 3; (5.30)

 

where A has been defined in (5.21) and

 

AT h
B = [8'9 ' (HT) (00" T+OTV0+BT)][Tm ‘ To ‘ (T '2']
M I
._AT AT-D__1._E Ell-
+[Bcg (Tm—“Go;+°cV0+BC)H’0M(TN 2+A 2 com A2 +M W10!
Vzgh
V0 AT AT h
-VO’+B- ATM ---—— 2-——(K+I<o)]
0 0 0 uOA (h )K: zfih W( ) 2 O T C 0
. _ ALT.
' [809 h ) (00Vp + OQVC + 80% DO (5.31)
2
a: v. +8 _ V09h(12)_1_(K MOHWW'.)
00 0 110(AT) hKO T c0 T COg
AT AT
(‘5‘) [ST + BCoO-oOwTwCoO) + vowT +ocoo)1 2}: (5.32)

82

 

 

.A9.
2 2gh
D=%- ZMAe B'g+ vooo+ 80 -—1——5-(AhT)(::AO:)
2a n sinh A? p
-(93)<ov+vo+e) (533)
h Co Op 0 °

Note that the zeroth—order solutions have also been used.
The solution to equation (5.30) which satisfies the

boundary condition (5.26) is easily found to be

 

 

_ g _ 2C 2 A2
(101 - A X:- + 32 + CZ + De (5.34)
%
h_h -g & +M
+ [D cosh (A7) -2-(B A +Ch sech(2 )]e (C100 aznh)

Finally the substitution of the temperature gradient
and the complete expression for 01 into equation (5.21) for

the composition gives

 

 

 

 

 

 

 

dgl = E + 2Fz + G doo + —§L— (5.35)
a wh
where
_AT _ Eh fl -51_h _
”(TI—{IT Tm To (h 2 + c("o h ( 2+A 2 C°th(A))
+fi-w -.. +0- m“ @552 M
M 10 p O O uO(AT) 11 0 aznh
+(93)2—p—(K +K o) (5 35)
h 2K0 T C 0 °
th
_1 AT2 _ o 1 .
F_§- (F) OT+ GCOO OO WR—g (KT +KCOO) (5.37)

 

 

83

 

 

 

 

 

 

 

 

 

 

 

 

(93 o - o - thO 53 (5 38)
h p 0 uOZAT) °
Equation (5.35) when integrated gives the follow—
ing expression for the composition:
C = R + Ez + Fz2 + E (5 39)
01 A '
where R is a constant of integration. The constant R is
found by substitution of (5.39) into the boundary condi-
tion (5.27) and using the previously found solutions for
C00, dOO' and do It is
h
E+o0 fiw) ): A2
we
R: h
2 sinh A?
h
T A—
F+° 0(773 B) (Ah-2)Ahe 2
' 2 2 ‘ h
A .2 sinh A5
' h
AT A—
+ °o(TT)C 6 _ (A2h2-3Ah+6)Ahe 2
A3 L_ 2 sinh A;
‘ h
AT h 3A—
Gong)D A2 h Ahe 2
+ e cosh A— - h
2A 2 Sinh AE-
- Mg coth A; (5.40)
2a n
where Q is
A92-
_ 2C h
Q — [D cosh(A% )— 2(3 — 7: + Ch) sech(A§)]e (5.41)

84

Equations first order in el.--The convective cor-

 

rections due to heat loss first enter in the terms of the

transport equations which are first order in 6. When the

terms of order 6 are equated, the following equations are

 

 

obtained:

v Vd + d + D4 V - v = 0

~10 ~ 00 00 2 ~ ~10

a nh
d + ”4 v - vc - v - J = o
00 2 ~10 00 ~ ~10
a nh
M

Ydlo + Ad105 + doo + ;§;; (8 YTlo + Voyclo)

 

 

: 2 s 1 .
BonoV Ylo + 80 (5“0 + ¢0)Y (Y Yio’

80(T +T00)A

 

 

 

 

 

 

 

2 M r AT
K v T = d + - c ——
o 10 ( 00 a2;; 36 vo(11
' 0 (AI) (Tr+TOO)B°V° - Q (A - v )
o h 85 0 ~ ~10
V gh B V
0 0 o
D0‘310 + ’TT’ ' “000 + ‘Eg‘ YTlo
v3 J v M -1
+ + u vc + d + va
3; o - 10 HT 00 22;; - 10
v —2
o M M
= —T d + a v a +
80 00 azflh 10~( oo azflh

 

 

 

(5.42)

(5.43)

(5.44)

(5.45)

(5.46)

 

1For the remainder of this section it will be assumed
that the quantities r and 2 have been scaled in such a way
that a and h may be taken to be unity. This greatly sim—
plifies the analysis to follow without any loss in gener-

ality.

85

It should be clear that some of the zeroth order solutions
have been used in these equations where convenient.

The boundary conditions for this order are

Tlo = 0 at z = 0 (5.47)
T10 = 0 at z = l (5.48)
BT10

W = - TOO at r = l (5.49)

for temperature,

 

 

2n 1 l
“f .f .[ dlordrdzdw = 0 (5.50)
O 0 0
2H 1 l
M _
‘[ .[ _[ [COOdlo + (100 + —§—— Clo]rdrdde-O (5.51)
a wh
O 0 0
for density and composition, and
vlOz = VlOr = J10z = 0 at z = 0,1 (5.52)
vlOz = vlOr = JlOr = 0 at r = l (5.53)

for the Velocity and diffusion flux.

Solutions first order in €.-—The increased number

 

and complexity of the equations for the convective correc-
tions necessitate a different approach from that used in
the previous two cases. In this case it is convenient to

begin with the solution for velocity. It is possible, by

86

introducing a velocity potential 4 which determines the

velocities according to

 

_ l M 8¢
Vior ‘ E (doc + ‘7" ‘2 (5°54)
a nh
and
-l
1 M ) a¢
V = - ... d + I (5.55)
102 r ( 00 azwh r

to satisfy the continuity equation (5.42) identically.
Although at first sight the equation of motion

(5.44) appears formidable, it can be used to obtain a

single partial differential equation for the potential ¢

alone. The curl of equation (5.44) has only one component:

 

 

 

 

 

 

 

 

 

 

A adlo + A d -+ M B aTlo + V aC10]
5r 00 ainh 0 8r 0 3r
3v 3v
2 lOr lOz
_ I _
‘ Bonev “‘52’ "—53- (5.56)
while the divergence of the equation of motion may be
written
2 a M 2 3
(v - A 55) le + d00 + V + A 53 (BOTlO + VOClO)
a nh
_ 4 2
‘ 86 (3“0 + ¢o) AV V102 (5'57)
. . 2 l 8
If equation (5.56) 18 operated on by V — —§ - Ag;

and equation (5.57) is Operated on by 5% , then each of the

87

left hand sides are equal and may be eliminated between
the two equations. This gives the following equation for

the velocities alone:

3v 8v

 

 

 

3v ¢
2 l. 3 2 lOr lOz _ 2 4 0 2_ l 102
V _ :3 A 5; V —§E__.._§E__._ A §‘+ FENV :7'—§E—— (5.58)

With the introduction of the velocity potential 4, equation
(5.58) becomes the following sixth-order partial differ-

ential equation:

 

-A2 2
2 l 3 2 l 3 .1 3 3 3
(V‘? Arab? ra—EEE”;;2'A52¢
-A2
2 4 ¢o 2 1 a J. a _
+ A 34+ “o))v - :7-e 55.; gr ¢ — o (5.59)

 

 

It is obvious that this equation is far too complex to be
amenable to exact analysis. There are however several
possibilities for approximations. If the terms of (5.59)
are expanded in detail, factors of the type (1 + A2f+ A4g+
...) occur very frequently. It was pointed out earlier
that A is typically on the order of 10-2. For this reason
the higher order terms in A are ignored wherever they occur.

This allows equation (5.59), after a great deal of algebraic

manipulation, to be written in the compact form

 

 

3 § L - (Eiin n3?¢ - + E § Ln' (”ilZA n35¢ ,==o (5.60)
n=l j=0 “3 r ar 3323 n=l j=0 3 r ar 3323
even odd

88

where the matrix Lnj is defined as

60 3 o o o o o
60 26 8 o o o o
30 22 15 4 o o 0

Ln. = (5.61)
3 10 12 11 12 1 o o
3 6o 6 12 3 6 o
1 o 3 o 3 o 1

Although equation (5.60) is still too complicated
to allow an exact analysis, an approximate solution will

be sought in the form

2 2 N M
¢* = r(r-l)z (z-l) r(r-l) + Z 2 C.
i=0 j=l 1

jzl sin wjr (5.62)

The approximating function has been chosen in such a way

that all the boundary conditions (5.53-53) are automatically
satisfied. Furthermore, each term is a member of a complete
set of functions and thus independent of every other term.
There are a great many techniques and criteria for the de-
termination of the coefficients Cik' but very little infor-
mation available which would allow a choice to be made from
among them. The method of weighted residuals (Finlayson

and Scriven, 1965 or Ames, 1965) has been successful in

solving a large number of complicated problems and was

selected for use in determining the Cik's.

89

The method of weighted residuals proceeds as follows:
Let the function R represent the residual obtained from the
substitution of the approximating function (5.62) into the

differential equation (5.60). The residual is given by

6 n n n 6 n n n
R = 2 2 Lnj (;EL :_?* 1 + 2 Lnj ('éinA anf: T (5.63)
n=l j=0 r Br 3323 n=l j=0 r Br 3323

 

 

even odd

The residual R is a measure of the extent to which the
function 4* satisfies the differential equation. As the
number of approximating functions is increased in succes-
sive approximations, it is hoped that the residual will
become smaller; the exact solution is obtained when the
residual is identically zero. As an approximation to this
ideal, the weighted integrals of the residual are set

equal to zero:
J) J) Wderdz = 0 j = 1,...NM (5.64)

If 4* is the exact solution, equation (5.64) is satisfied
regardless of the choice of weighting functions.

The weighting functions,W3, can be chosen in several
different ways, and each choice corresponds to a different
criterion in the method of weighted residuals. The least
squares method, which apparently was first presented for

this type of application by Picone (1928), uses the

90

weighting functions 3%3— . The corresponding interpretation
ik

is that the mean square residual

=f [derdz (5.65)

is minimized with respect to the constants Cik‘ (See
Mikhlin, 1964), for a discussion of the least squares
method.)

In practice the Cik's are determined by a direct
minimization of the integral (5.65) using a computer rou-
tine developed in connection with least squares statistical
problems. Of course the numerical values of the Cik's
depend on the Specific problem, in particular on the value
of A. In Appendix C, numerical values calculated for a
typical value of A are presented together with plots of
the velocity profiles resulting from equations (5.54-55).

Inspection of the temperature equation (5.45) re-
veals that once the velocity is determined, the only
dependent variable in the equation is the temperature
correction T10. Introduction of the velocity potential
from equation (5.55) yields the following equation for

the temperature corrections:

2
aT10 + a T10 1 84*

r 322 = (plz + p2) FEB—r— (5.66)

 

 

_3._r
8r

HIP

9l

where
80(A'I‘)
p1 = 13.8.5— A - oovomT) (5.67)
and
= i. (o o - c )(AT) + 8° T - l<AT)) (A—o v (AT) (5 68)
92 K0 0 0 v0 63' m 2 o 0 °

The 2 dependence of T10 is most easily determined by taking
a finite sine transform of equation (5.66) and the asso-
ciated boundary conditions (5.47-49). Multiplication of

(5.66) by sin(nnz) and integration over 2 from 0 to 1 yields

%»5% r d? - n NZJT (r) = Thn(r) (5.69)

after two integrations by parts and use of the boundary

conditions. The sine transform of T10, denoted by T in

n
(5.69) is defined as
l

Tn(r) if Tlo(r,z) Sln nnzdz (5.70)
0

while the inhomogeneous term Thn is given by the following

integral:

0)

_ l 4* .
Thn(r) e]. (plz + p2) E 3r Sln nnzdz (5.71)

 

92

The boundary condition (5.49) is transformed into

8T
n _’ 1 _ _ _ n _ n
.5? _ H? (Tr Tm + 1/2 AT)(1 ( 1) ) + ( l) (5.72)

Finally, because of the completeness and othogonality of

the sine functions on the interval 0 to 1, T is deter-

10
mined from the transforms Tn, by the sine series

I (r) sin (nnz) (5.73)

(r,z) = 1 n

n

HMS

T10

The homogeneous part of equation (5.69) is imme-
diately seen to be a special case of the modified Bessel

equation. The general solution can be written as

T (r) = ClIo (nnr) + C (nwr) (5.74)

n 1K0

where I0 and K0 are the modified Bessel functions of order
zero, of the first and second types, respectively. The
function K0 (r) can be approximated by -ln (r) for small r,
and thus increases without limit as r approaches zero.
Consequently the coefficient C2 will be set to zero. A
particular solution to equation (5.69) is much more diffi-
cult to obtain. One procedure which works well in an
approximate sense is to write the inhomogeneous term as a
power series in r which is easily derived from the expan-
sion of the sine functions in 4*. If this series is rep-

resented by

where

and

The first few coefficients of (5.75) are

n0

n1

n2

where

T = Z

“P 2:0
Gno = Gnl '
(2 + 2)2 G

N M
Z ,1
i=0 j=l
N M
-,Z .2
1=0 3=1

fljC

n,£+2

ij (lei+1,n

93

- nznzG = H

- 2 (p1K1,n + p2K0,n)

the matrix Km n is defined as

X

j=0

I

(-l)j (2j)l(mN)-2j-l[

n,£ n,£

n+4
Zj

(5.75)

(5.76)

(5.77)

(5.78)

(5.79)

+ PZKin)-6(p1Kl’n+p2Ko’n) (5.80)

ij(lei+1,n+p2Kin)+4(lel,n+pZK0,n) (5°81)

)-2(3§3J+(3§2))

94

The coefficient C is to be determined from the

1
boundary condition (5.72) and has the following value:

c = 1 (Tr-Tm + 1/2 AT)(l-(-l)n) (5.83)

l (nn)2Il(nw)

 

+ (- 1)n - nn 2 26

2 n2

The temperature correction has thus been com-
pletely specified and may be calculated for any particular
set of conditions. In Appendix C the temperature profiles
have been calculated and graphed for the same conditions
under which the velocity corrections were calculated.

If the potential V for the diffusion flux is de-

fined according to

£10 = 34 (5.84)

then the equation of continuity (5.43) becomes

34*

2 l
0 (AT) '1: '51.— (5.85)

V W = - o

This equation is virtually identical to equation (5.66)

for the temperature. The solution to (5.85) is therefore

W = 2 f (r) sin (nnz) (5.86)

where

(D

f (r) = B I (nnr) + Z GAR:Q (5.87)

2:2

95

where the equations (5.77-82) are again used to determine

Gig but with pl equal to zero and p2 = - 0 (AT). The

0
constant B1 is determined from the boundary condition
(5.53) and is equal to

_ -1
Bl - EETHFT :2 Gn2 (5.88)

With the velocity, temperature, and diffusion flux
corrections completely determined, it is now possible to
obtain solutions for the composition and density correc—
tions from the diffusion flux equation (5.46) and the equa-
tion of motion (5.44). It is a simple matter to rearrange
and eliminate terms between these two equations to yield

the following two equations:

 

 

 

 

 

 

VODO vggh )
(Y + A5’d1o ’ doo + ‘2}; ‘6" ~10 + 636T ' v00 0 B0 ”1%.
+(g B' W)[ 2vl + (1/3 no + ¢o)Y(Y'Y1o) (5.89)
Vogh D0
Yc1o = 0o ‘ FEET YT10 ‘ E; 910
v (d M )‘1
_ o oo znh 2 .
“0 [V Ylo + (1/3 ”0 + ¢O)Y(Y y10% (5.90)

The right hand sides of both of these equations are now
completely known and it is an easy matter to integrate

the equations. However, the major interest in the thermal

96

diffusion experiments usually centers on the gradients,
and it is therefore deemed unnecessary to go further.

This completes the solution to the convective heat
loss equations at first order. It is obvious that the same
process could be carried on to higher orders, but inasmuch
as the corrections at the first order are already quite
small, the corrections from higher orders are expected to

be negligibly small.

6. Summary

 

It has been shown in this chapter that the pertur-
bation expansion method can be a useful and valid technique
for the solution of complex partial differential equations
under a variety of difficult conditions. Corrections to
account for convection due to heat loss in a thermal dif—
fusion cell are now available for the first time, as well
as corrections for the variability of phenomenological
coefficients. It is doubtful that another technique could
have achieved the same results without much greater effort.
Finally, because of the simplicity of the approach, it
seems reasonable to expect that it can be applied easily

to countless other difficult problems.

CHAPTER IV

PARAMETER ESTIMATION

1. Introduction

 

The preceding chapters have been concerned with the
problem of determining theoretical descriptions for experi-
ments on liquid mixtures. The approach has been as rig-
orous as possible even though practical necessities have
at times required approximation methods such as the per-
turbation technique. Once the differential equations
appropriate to a particular experiment have been formue'
lated and solved, the next step is to compare the solutions
with actual experimental data.

Much of modern experimental physical science is
concerned with the assignment of numerical values to quan-
tities which occur in physical theories. This assignment
is never absolutely accurate. No matter how well a mea-
suring device is designed, there will always be limita— -
tions inherent in the materials, methbds of construction,
or human observation which will make the reproducibility
of a measurement impossible at some level.

The inability to make precise measurements is only

one of the problems which arises when the results of

97

98

experiments are to be tested against physical theories.
Another problem which is becoming more and more important
as the sophistication of experiments and theories increases,
arises in the fact that the quantities which can be most
conveniently measured are not the quantities of basic
interest in the theory. Mathematical relationships must
be used to go from one to the other; these relationships
may be simple or complex. In any case, if the quantity

A has been measured with a given precision, we will want

to know how well we know the value of a derived quantity

B. A more general problem arises when a multitude of quan-

tities A1, A ..., An have been measured, and it is nec-

2:
essary to determine the values of a set of theoretical
parameters 31' 32, ..., Bm all of which are interrelated,
each being dependent on all the Ai'

How are the best values of the Bi to be determined
and is it possible to make a reliable estimate of the un-
certainty in each estimated parameter? In fact is it
possible that some of the parameters are not determined
at all from the data? The solutions to these and related
problems are the concern of statistical analysis. The com-
parison of experimental and theoretical descriptions can
only be made after an attempt to answer some of these
questions by statistical analysis.

In cases for which the relations between the experi-

mental data and the parameters to be estimated are linear,

99

standard techniques of multiple regression analysis have
been developed to answer some of the questions above. If
the random errors involved are normally distributed, then
many of the standard statistical tests (F test, Student's

t test, 222') provide very powerful methods for testing
hypotheses on the parameters. Unfortunately, the nonlinear
problem occurs quite frequently, and typically the distri-
bution of errors is unknown. In these cases statistical
analysis has provided very few methods of general appli-
cability.

In this chapter a generalization of the least
squares or minimum variance parameter estimation technique
is presented (see Hamilton, 1964, or Wentworth, 1965).

The present approach is novel in several respects: (1)

No assumptions are made about the data being uncorrelated;
(2) It is not necessary to assume that the eXperimental
errors are normally distributed; (3) No assumption of
linearity in either the data or the parameters is necessary.

Finally, two examples which arose in the study of
pure thermal diffusion experiments and which illustrate
the methods and the computer techniques necessary to imple-
ment them are presented. Of particular interest in these
examples is the difference between results obtained from
the new approach and those obtained by the ordinary least

squares method.

100

2. Generalized Least
Squares Adjustment of Data

 

 

It will be supposed that in a given experiment or set

set of experiments the n quantities

Z 22’ 00., Zn (201)

1'

are measured and further that each of these n observables

is known to depend on a set of m parameters
Y2, ..., Y (2.2)
as represented by the function

f(Zi, ..., Z*

n, Y* .00, Y*) = 0 o (2.3)

l' m

The symbol * in (2.3) denotes the true or theoretical
value. The value of the ith observable for the kth ob-
servation will be denoted by zik' Thus the data for a
series of N experiments or N repetitions of an experiment

can be represented as

Z Z ooo'Z - I i=1, ...,N (204)

li’ 2i’ n1

It will be assumed that the errors or deviations of
the experimental measurements from their theoretical values
are random variables having a joint distribution with zero

means, and that the variance-covariance matrices

101

-1

(Ml£)kj - z* (2.5)1

Z

_ *
(Zki j2 ij

are of rank n. Thus one of the basic assumptions of this
approach is that the errors in the observations are drawn
from p0pulations with finite second moments. It may be
noted at this point that this is the only requirement,
and that the theory of least squares does not depend on
the assumption that the errors follow a normal distribu-
tion, although the latter assumption is commonly and mis-
takenly believed to be necessary.

A good measure of precision is usually taken to
be inversely proportional to the dispersion or variance
of a random variable. Thus one method of adjusting the
data might be to choose the estimates of the theoretical
data points, Zij’ in such a way that the total sample

variance

 

1Consider a sample space, a random variable x de-
fined over this space, and some function of x: f(x).
The expected value of f(x) is defined as

E {f(x)} = 2 f(x) P(x)

if the sample space is discrete, with probability measure
P(x): the sum is over all points of the sample space. If
the sample space is continuous with probability density
function ¢(x), the expected value 1A; defined as

E {f(X)} = fmf(x)4(X) dx

—oo

In either case, the expected value is defined only if the
sum or integral exists.

102

N N n n ..
s = Z Z 2 Z (zki -Zfii) Mia (sz - ZEj) (2.6)

i=1 j=l k=l =1
is a minimum. Note that in calculating the total sample
variance the experimental points have been weighted in
inverse proportion to their variances. Unfortunately, the
estimates Zij cannot all be varied independently in the
minimization of (2.6) since they are related to each other

and to the parameters by the equations

f (

*
i zli'z* zai’Yi'oooY;l) - 0 o (2.7)

21,000

It is usually convenient, when minimizing under a con-
straint like (2.7) to introduce Lagrangian multipliers,

A. and to minimize

1'

N
I ...
s — s + E 11f. (2.8)

instead of (2.6). Equating to zero the Ar' Yfi, and 2;:

derivatives of 8' yields the minimum conditions

 

 

 

fr = O r=lloooN (2.9)
N air

2 Ar 5?? = 0 k=lpooom (2010)
r=l k

and
8f. 31 .2. r.
l A T) = M (Z . 2*.
2 r azpr i=1 k=1 pk k1 k1
r=l, .N

p=l,...n (2.11)

103

Note that the symmetry of Mik has been used in simplifying
equation (2.11).

Equations (2.9-11) are N(n + l) + m conditions for
the determination of the N(n + l) + m unknowns Ar’ Yfi, and
Zij‘ Unless the equations of (2.9) are linear, these are
the most difficult to solve. In the nonlinear case the
standard technique of Newton's method may be used, result-
ing in an iterative procedure. Newton's method for (2.9)
consists of expanding each equation in a Taylor series
about some initial guess for the parameters and truncating

after first derivatives:

 

 

 

 

o .2. 3f. o If 3f.
f = f " (Y. " Y5) " ——T— (Z. "Z? ) = 0 (2.12)
r r j=l SY; j j j=l azjr jr jr

If the differences (Yg - Y3) obtained from the solution of
(2.10-12) are sufficiently small, then the Yg's may be
taken as the least squares estimates of the parameters Y3.
If the differences are not small enough, then new estimates
of the parameters can be found by subtracting the differ-
ences from the original estimates and re-solving the equa-
tions. This procedure may be continued until the differ-
ences are as small as desired.

Terms containing the differences (Z.

_ * '
jr Z jr) Wlll

be eliminated from equation (2.12) by using the solution

to equation (2.11). Equation (2.11) yields

104

 

 

N n ._
(zjr ' Zgr) = 5 2 Hri A ifé—- J_l'...n
1:1 k=l 3k 1 azki r=l,...N

where Hg: is determined from the conditions
1:]. N MR, Hri - 6 6 21:1,000N

._ mj jk ‘ 12 km _
3"]- Fl m-lpooon

If equation (2.13) is introduced into (2.12) one obtains

0 m 31:r 0 1
fr - jgl Ff? (yj - 2;) — 7 Z (w )rixi (2.15)

where for convenience the following definition has been
made:

3f.
1

azki

Bfr Hri i=1,...N
Ber jk

 

 

 

 

r=l,...N

Equation (2.15) may be written in matrix form as follows:

 

 

F0 - AAY = % w’lx (2.17)
where
o _ o
(F )k — fk (2.18)
(AY) = (Y0 - Y*)
k k k
Bfi
(A)ij = BY;

(2.13)

(2.14)

(2.16)

(2.19)

(2.20)

105

and the matrix W has been defined in (2.16). Equation

(2.10) in matrix notation becomes

ATA = o (2.21)

Multiplication of (2.17) by ATW and use of (2.21) yields

AY = (ATWA)-l AiWFO (2.22)

which corresponds to the normal equations usually obtained
in the ordinary least squares method. Thus to obtain
estimates of the deviations of the parameter guesses from
their true values, it is only necessary to perform matrix
inversions and multiplications.

Suppose now that (Mij);i is known only to within a

scale factor:

ij -1 _ 2 ij -1
(M )k2 - o (N ’k2 (2.23)

It therefore follows that W is determined to within a

constant:

C

where P is determined from (2.14) and (2.16). Using (2.24)

in (2.22) gives

AY = (.zx.TP.zx.)‘l ATPFO (2.25)

~ ~

and the normal equations are thus seen to be independent

of the scale of the weight matrix P.

106

Errors in the parameter estimates.--Once the best

 

least squares values have been obtained as solutions to
(2.25), it is necessary to determine the variance-covariance

or moment matrix for the parameters. The desired matrix is
Q =e{(§Y)(gY)T} (2.26)

By use of the normal equations (2.25), (2.26) may be

written

0 = (ATWA)-l ATWE {gogoT}WA(ATWA)'1 (2.27)
Now if the iterative procedure has converged so that the
deviations are small then

s {ggT} = w“1 (2.28)
Use of (2.28) in (2.27) gives

Q = (ATWA)-l = 02 (ATPA)-l (2.29)

To obatin the moment matrix for the parameters, it is thus
necessary either to know 02, in which case Q is completely
determined, or to make an unbiased estimate for 02, in
which case only an unbiased estimate of the moment matrix
is available. It will now be shown that an unbiased esti-
mate of o2 is indeed available from the results of the

least squares treatment.

107

The value of S in equation (2.6) at the least
squares solution is found by substitution of (2.13) into

(2.6), with the result
s = (F0 - AAY)T w (F0 - AAY) . (2.30)

Use of the normal equations (2.22) in (2.30) gives

T
s = g0 w§° - (gY)T(ATWA)(§Y) (2.31)

Because of (2.28-9) the expected value of (2.31) is
e{S} = N - m (2.32)

From (2.32) and (2.24) it is clear that an unbiased esti-

mate of 02 may be obtained from

2

o = (N - In)—1

(F0 - AAY)TP(FO - AAY) (2.33)

which means finally that an unbaised estimate of the moment

matrix of the parameters may be obtained from

Q = (N - m)-l (F0 - A9Y)TP(§° - AéY)(ATPA)—1 (2.34)

3. An Alternative Approach

 

It happens very often that the iterative procedure
suggested in the previous section fails to converge because
the initial parameter estimates are not good enough. This
is typical behavior for analogues of Newton's method.

Several procedures have been suggested for improving the

108

convergence of the least squares method. Box (1958),
Hartley (1961), and Strand, Kohl, and Bonham (1963) have
all suggested the introduction of a parameter t into the

normal equations

AY = t (ATWA)-l ATWFO (3.1)

~

A value for t is then chosen which makes the variance S

of equation (2.6) as small as possible. While this pro-
cedure extends the range of convergence of the method,
divergence is still obtained for a large number of prac-
tical problems, which necessitates even further refinement
of the method.

The method of damped least squares as suggested by
Levenberg (1944) and subsequently modified by Meiron (1965)
has been one of the most successful improvements in the
basic least squares method from the standpoint of conver-
gence. (See the comparisons of Pitha and Jones (1966).)
This method requires that the changes of the variables
éYTGéY be minimized simultaneously with S. This keeps
the method from taking steps outside of the region of
validity of the Taylor series expansions, one of the more
important reasons for the divergence of the ordinary least
squares method. Meiron's particular choice of G was to
take the diagonal elements of (ATWA). The parameters are

then determined from

109

AY = (A¢WA + t-.G)'1 A¢WFO (3.2)

~

where again the parameter t is to be chosen so as to make

S as small as possible. While the damped least squares
method extends the range of application of the least squares
technique, it suffers from two computational disadvantages:
(l) The derivative matrices must be calculated at each step
of the procedure. (2) A large number of matrix inversions
are required not only in determining the parameter esti-
mates, but also in determining the optimum values of the
parameter t. These disadvantages can often increase sig—
nificantly the amount of time required for convergence to
be obtained.

For these reasons an alternative approach was
sought which would eliminate some of the convergence and
computational difficulties encountered in the ordinary or
damped least squares methods. An approach which was found
to be successful is summarized as follows: The minimiza-

tion of the quantity

8" = F'T w F' (3.3)
with respect to Yfi results in the following minimum condi-
tions

3F. T I _

52% w E 0 , (3.4)

 

110

where W is the matrix defined in (2.16) and F' is obtained

from the Taylor series

I
n 8f.
' * * = 1 - *
fi(zli'22i,...zni’Yl'...Ym) jgl F2; (Zji Zji)+ .0. o (3.5)

If the behavior of the function F' near the minimum can be

represented accurately by the approximation

*3. = 50 + 1‘9“ (3.6)

then the values of the parameters Yfi which give a minimum
for S" are the same values as those predicted by the nor-
mal equations of the least squares method.

The minimum conditions (3.4) are not used in this
case for the determination of the minimum. Instead a one
dimensional searching technique, which operates directly
on the function (3.3) and which produces a sequence of
conjugate directions for the searches as it proceeds, is
used. This method, first introduced by Powell (1964),
requires no derivative evaluation or matrix inversion.
Furthermore, because of the use of conjugate directions
the method can find the minimum of any quadratic function
of n variables in exactly n searches. The algorithm given
by Powell has been programed and the listing of the pro-
gram and directions for its use are given in Appendix B.

The next two sections illustrate the use of the

method on two nonlinear functions of importance in

111

analyzing pure thermal diffusion experiments. It is impor-
tant to note the effects of different choices for the

weight matrix W.

4. Fringe Shape Analysis

 

In his experimental determination of thermal dif-
fusion constants for the CCl4 - C6H12 system using wavefront-
shearing interferometery, Anderson (1968) observed inter-
ference fringes of skewed parabolic shape. The curvature
of these fringes is due principally to the temperature
dependence of the index of refraction of the CCl4 - C6H12
mixture. This curvature is also found when a temperature
gradient is applied to either pure component. The problem
then is one of relating the curvature of the fringes to
the temperature dependence of the index of refraction.

The curvature of the fringes is related to the

index of refraction by the following formula:
x=A (n(2+bl) -n(z -bl)) +8 (4.1)

where A and B are constants, n is the index of refraction,
x is the displacement of the fringe measured from an
arbitrary vertical plane, 2 corresponds to a vertical
distance in the thermal diffusion cell, and b1 is a shear-

ing distance introduced by the interferometer.1

 

1The formula derived by Bryngdahl (1963) and sub-
sequently used by Anderson (1968) has the same appearance
as (4.1). However, the constant A, as they define it,
makes their equation dimensionally inconsistent. The cor-
rect definition is A = a- df/l , where a is the path length
through the cell, df is the distance between fringes, and
A is the wavelength of light used.

112

The temperature distribution in the cell is almost
always well known. For example in the case of a pure com-

ponent, the temperature is given by

ATz + Tm - % KT (AT)2 (22 - 1) (4.2)

a
u
NH“

when a temperature difference AT is applied between the
horizontal planes bounded by z = — l and z = 1. The
averageof the temperatures on the two planes is Tm, and
KT denotes the temperature variation of thermal conduc-
tivity, K:

O)

K

'55) (4.3)
It has also been assumed that the thermal conductivity
depends only linearly on temperature and that the tempera-
ture variation is entirely vertical. The index of re-

fraction will be represented by the following series:

 

 

 

 

E T-Tm k
n (z) = u —— (4.4)
The derivatives of index of refraction with respect to
temperature are determined from
i 'l
L3). = .317 ”i (45,
ST (AT)
T
m

once the coefficients ui are known.

113

Substitution of the index of refraction-temperature
expansion (4.4) into the interferometer equation (4.1)
gives the final working equation for the least squares

analysis:

fi = xi - Ago uk(AT)‘k (T(zi+bl)-Tm)k - (T(zi-bl)-Tm)k - B (4.6)
where the fi has been written to correspond to the notation
of the previous sections. It will be assumed that the
data, xi and 21, are uncorrelated and of equal weight for
all i. Although experience with the measuring technique
is limited, nothing has been observed to date which would
invalidate this assumption.

The coefficients are determined according to the
least squares analysis from the minimization of the quan-
tity

N

Z f.W. .f. = F WP (4.7)
1 j=1 1 ij j 2 ~

0
II
II M2

i

It is particularly important to determine the weight
matrix W carefully according to equation (2.16). One
finds all too often that in problems where the data are
uncorrelated and of equal weight, the weight matrix has
been taken to be the identity matrix. This can result in
serious systematic errors in the estimates of the coef-
ficients for some problems; the fringe analysis problem
provides a good illlustration of errors that can result

from improper weighting.

114

Analysis of the fringe shapes has been done auto-
matically using the computer routines FRNGFIT and FRNGSTAT
in conjunction with the MINIMIZE routine. These programs
are listed in Appendix B together with instructions for
their use. Analysis of the fringes observed for the C6H12

system maintained at a mean temperature of 24.2°C with a

3.382°C temperature difference gives the estimates

_ _ -5
u2 — 1.9 x 10
u3 = 8.1 x 10"5

_5 (4.8)
114 = 7.1 x 10
u5 = 5.2 x 10"5

under the assumption of unit weights, while use of weights

calculated from (2.16) gives the following estimates

u2 = - 2.2 x 10‘5
u3 = 9.0 x 10'5
-5 (4.9)
u4 = 3.5 x 10
u5 = - 1.5 x 10‘5 .

Estimates of the standard deviations of these coefficients

are
02 = 1.1 x 10'6
03 = 3.6 x 10-6
-5 (4.10)
04 = 1.1 x 10
o = 2.6 x 10‘5

115

The first coefficient ul is determined in a separate
experiment and is therefore not estimated here.

It is apparent from the magnitudes of the standard
deviations that only the first two coefficients are deter-
mined with satisfactory precision from the fringe shapes.
Even for these two estimates, however, the difference be-
tween the estimates in the two cases is on the order of
three or four standard deviations, which clearly demon—

strates the need for proper weighting.

5. Time Constant Analysis

The preceding example illustrated the need for
careful choice of the weight matrix. However, the fringe
shape equation (4.6) is 1iinear in the coefficients and
thus does not illustrate the advantage of the MINIMIZE
procedure over the analogues of Newton's method. In this
section another least squares problem from the analysis
of thermal diffusion experiments is presented. In this
case the working equations are not linear in the parameters
and therefore provide a good test of the MINIMIZE routine.

Anderson (1968) observed that during the course
of a thermal diffusion experiment the fringes moved hori-
zontally with time. The time dependence is exponential

and the fringe shift with time can be expressed by

d = ae—bt + c (5.1)

116

where d is the fringe shift measured from an arbitrary
initial reference plane, t is time, and a, b, and c are
constants to be determined by least squares analysis.
Anderson showed that the ordinary diffusion con-
stant could be calculated from the value of b while the
thermal diffusion factor could be calculated from the value
of a. The computer routines EXPO and EXPOSTAT listed in
Appendix B were written for the purpose of estimating the
constants a, b, and c from the experimentally measured

fringe shift with time according to the equation

f. = d. - ae_bti - c . (5.2)

To illustrate the power of the method, the fringe
shift data reported by Anderson (1968, p. 158, Table Sj)
was analyzed with these routines. The minimum was easily
found to within 0.01 percent in only four or five itera-
tions, using at most four seconds of computer time. The

estimates found are

a = - 7.11
b = 1.81 x 10"2 (5.3)
c = 1.01 x 101

and the estimated standard deviations are

0 = 1.5 x 10‘1
a

db = 5.6 x 10'4

o = 3.3 x 10‘2

117

The smallness of the standard deviations indicates that
the analysis of the time shift behavior can produce re-
liable estimates of both the ordinary diffusion constant

and the thermal diffusion factor.

6. Summary

 

It is apparent from the results of the last two
sections of this chapter that the problem of estimating
the values of interrelated theoretical parameters from
correlated experimental data is within the means of most
physical scientists. The general theory of parameter
estimation set forth in Section 2 has been reduced to
practice and made efficient by the use of a novel minimi-
zation scheme. The method can provide estimates for a
wide range of complicated problems and, simultaneously
gives estimates of standard deviations and goodness of

fit information.

CHAPTER V
CONCLUSION

In the preceding chapters a general approach to
the principles and practices of nonequilibrium thermo-
dynamics of fluid mixtures has been presented. Examina-
tion of previously existing theories of irreversible
processes revealed many inconsistencies, especially with
regard to phenomenology and convective behavior. The
reformulation of the theory given in Chapter II embodies
a new approach to the specification of phenomenological
relations through the rational use of the entropy in-
equality. One of the most significant improvements in
this respect is the prOper enumeration of independent
variables. The approach also resolves such questions as
(l) validity of a Gibbsian differential equation in non-
equilibrium situations, (2) neglect of inertial and
viscous terms in the phenomenological equations, (3)
neglect of the kinetic energy of diffusion, (4) defini-
tion of mechanical equilibrium, and (5) appropriateness
of the Newtonian stress tensor for mixtures.

A particularly interesting and novel result is
the requirement that the phenomenological equations for

diffusion fluxes are determined by componentwise equations

118

119

of motion. The latter contain the partial stress tensors,
which are phenomenologically related to the componentwise
velocity gradients by as many as ten viscosity coefficients
for a binary mixture. Snell and Spangler (1967) have also
deduced a more general relation between stress and velocity
gradients for mixtures. Although their macroscopic treat-
ment is less extensive and uses friction coefficients
rather than viscosity coefficients, the conclusion that a
more general relation than the Newtonian stress formula is
appropriate for mixtures is substantially the same as that
reached here. This is a promising area for new experimental
study. Hopefully, the viscosity coefficients (or the re-
laxation times to which they are related) can be deter-
mined from light scattering or ultrasonic measurements on
mixtures under the influence of temperature and concen-
tration gradients. Moreover, there are interesting

aspects of microscopic theory (Bearman and Kirkwood,

1958; Snell, Aranow, and Spangler, 1967) to reconsider

now that there is a firm macrosc0pic foundation for the
existence of partial stress-strain relations.

In order to make use of the partial differential
transport equations derived in Chapter II it was necessary
to develop better mathematical techniques for their solu-
tion. Chapter III introduces the simple but powerful
pertubation expansion method and illustrates its use in

the solution to convection and variable phenomenological

120

coefficient problems of pure thermal diffusion-experiments.
Even though this application of the method is fairly spe-
cific for illustrative purposes, it should be clear that
the method is general enough to be applied to a vast number
of other problems.

Besides the perturbation method, the method of
weighted residuals (in the least squares form) and the
finite transform method were used effectively in the
determination of the velocity and temperature profiles
for the heat loss problem of pure thermal diffusion. The
use of these three methods in tandem provides an approach
so powerful that the solutions to any complicated system
of partial differential equations are now within reach.

It should be mentioned that the approximating
functions, expansion parameters, and transform functions
used in Chapter III are in some sense arbitrary. It re-
quires a great deal of experience and cleverness to make
the choices which give Optimum results. For this reason
the results should not be thought of as final, but sub-
ject to refinement by a better choice of procedures.
Also, in order to give practical answers, some of the
expansions used in Chapter III were truncated without
investigation of the convergence of the series, perhaps
resulting in inaccuracies. These areas for improvement
are mentioned only in the hope that future application of

these methods will be even more successful.

121

Finally the important subject of parameter estima-
tion and comparison of theory with experiment has been
treated in Chapter IV. For the first time a detailed
account of the theory of estimation of nonlinearly inter-
acting parameters from data with correlated eXperimental
errors has been presented. The reduction of the estimation
theory to practice has been accomplished with the intro-
duction of a minimization scheme which is computationally
more convenient than the usual least squares method with
its associated matrix inversion problems. The illustra-
tive examples from pure thermal diffusion experiments do
not test the full power of the method, and it is certain
that much more difficult problems can be easily solved.
Convergence of the computational algorithm is probably the
most difficult remaining problem in the parameter estima-
tion field and any contribution to the state of the art
in this matter will be extremely valuable. The interest-
ing new approach of Kizner (1964) may ultimately lead to

improvements in the method described here.

BIBLIOGRAPHY

BIBLIOGRAPHY

Agar, J. N. and J. C. R. Turner. Proc. Roy. Soc. (London)
225A, 307 (1960).

Ames, W. F. Nonlinear Partial Differential Equations in
Engineering (Academic Press, New York, N.Y., I965).

 

 

Anderson, T. G. Ph.D. Thesis, Michigan State University,
East Lansing, Michigan, 1968.

Bearman, R. J. and J. G. Kirkwood. J. Chem. Phys. 28, 136
(1958). "

Bellman, R. Perturbation Techniques in Mathematics, Physics,
and Engineering (Holt, Rinehart, and Winston, Inc.,
New York, N.Y., 1964).

 

 

Bowen, R. M. Arch. Rational Mech. Anal. 24, 370 (1967).

Box, G. E. P. Bull. Inst. Intern. Statistique 36, 215
(1958).

Bryngdahl, o. J. Opt. Soc. Am. 53, 571 (1963).

Coleman, B. D. and W. Noll. Arch. Rational Mech. Anal.
13, 167 (1963).

and V. J. Mizel. Arch Rational Mech. Anal.
13, 245 (1963).

and V. J. Mizel, J. Chem. Phys. 40, 1116 (1964).

Curie, P. Journal de Physique, 3d series, 3, 393 (1894.

de Groot, S. R. L'Effect Soret (North Holland Pub. Co.,
Amsterdam, 1945).

 

. Thermodynamics of Irreversible Processes (North
Holland Pub. Co., Amsterdam, 1952).

 

and P. Mazur. Nonequilibrium Thermodynamics
(North Holland Pub. Co., Amsterdam, 1962).

 

122

123

Finlayson, B. A. and L. E. Scriven. Chem. Eng. Sci. 33,
395 (1965).

Fitts,ILlD. Nonquilibrium Thermodynamics (McGraw-Hill,
New York, N.Y., 1962).

Green, A. E. and P. M. Naghdi. Arch. Rational Mech. Anal.
33, 243 (1963).

. Int. J. Eng. Sci. 3, 231 (1965).

Haase, R. Thermodynamick der Irreversiblen Prozesse
(Stéinkopf Verlag, Darmstadt, Germany, 1963).

Hamilton, W. C. Statistics in Physical Science (Ronald
Press, New York, N.Y., 1964).

Hartley, H. O. Technometrics 3, 269 (1961).
Horne, F. H. J. Chem. Phys. 33, 3069 (1966).
and R. J. Bearman. J. Phys. Chem. 33, 4128 (1967).

Hurle, D. T. J., J. B. Mullin, and E. R. Pike. J. Chem.
Phys. 13, 1651 (1965).

Ingram, J. D. and A. C. Eringen. Int. J. Eng. Sci. 3,
289 (1967).

Kirkwood, J. G. and B. Crawford, Jr. J. Phys. Chem. 33,
1048 (1952).

Kizner, W. J. Soc. Ind. Appl. Math. 33, 424 (1964).

Levenberg, K. Quart. Appl. Math. 3, 164 (1944).

Meiron, J. J. Opt. Soc. Am. 33, 1105 (1965).

Meixner, J. and H. G. Reik. Handbuch der Physik, Vol.
iggéf, S. Flugge, Ed. (Springer-Verlag, Berlin,

Mikhlin, S. G. Variational Methods in Mathematical Physics
(Pergamon Press-MacMillan, New York, 1964).

Mfiller, I. Arch. Rational Mech. Anal. 33, l (1968).
. Arch. Rational Mech. Anal. 26, 118 (1967).

Noll, W. Arch. Rational Mech. Anal. 3, 197 (1958).

124

Onsager, L. Phys. Rev. 33, 405, 33, 2265 (1931).

Picone, M. Rend. Circ. mat. Palermo 33, 225 (1928).

Pitha, J. and R. N. Jones. Cand. J. Chem. 33, 3031 (1966).
Powell, M. J. D. The Computer Journal 3, 155 (1965).

Prigogine, I.‘ Etude Thermodynamique des Phénoménes Irre-
versible (Editions Desoer, Liége, 1947).

, Introduction to Thermodynamics of Irreversible
Processes (C.C:Thomas, Springfield, Ill, 1955).

. Physica 33, 272 (1949).
Smith, G. F. Arch. Rational Mech. Anal. 33, 182 (1965).

Snell, F. M. and R. A. Spangler. J. Phys. Chem. 13, 2503
(1967).

, R. Aranow, and R. A. Spangler. J. Chem. Phys.
33, 4959 (1967).

Strand, T. G., D. A. Kohl, and R. A. Bonham. J. Chem.
Phys. 33, 1307 (1963).

Thomaes, G. Physica 33, 885 (1951).
Truesdell, C. Rend. Lincei (8) 33, 33 (1957).

and R. A. Toupin, Handbuch der Ph sik, Vol. III/1,
S. Flugge, Ed. (Springer-VerIag, BerIin, 1960).

and W. Noll, Handbuch der Physik, Vol. III/3,
S. Flugge, Ed. (Springer-Verlag, Berlin, 1965).

Wentworth, W. E. J. Chem. Ed. 33, 96 (1965).

APPENDICES

APPENDIX A

RESTRICTIONS ON PHENOMENOLOGICAL COEFFICIENTS

IMPOSED BY POSITIVE ENTROPY PRODUCTION

It was suggested in Chapter II that the requirements
of non-negativity of entropy production imposed by the sec-
ond law of thermodynamics not only restricted the magnitudes
of the phenomenological coefficients, but in addition served
to limit the number of independent coefficients.

It is usual in the theory of nonequilibrium thermo—

dynamics to write the entropy production in a bilinear form:

11'
o=)g-xzo (Al)

and to write phenomenological relations for the fluxes Ja’

and the forces Xa’ in the following way:

a=0,...fl , (A2)

a
9
II
II M=l

Q
a

x
0 B”

B B

where the (0+1)2 quantities 9&8 are the phenomenological
coefficients. Often it is possible to show that a relation

of the type
‘2’

Q = O , 8:0,...0 , (A3)
0:1 as

125

 

126

exists among the coefficients. This implies that only

n (n+1) of the original (n+1)2 coefficients are independent.
One important example of this occurs when the fluxes JG,
a=1,...v, are the diffusion fluxes, related by

V

J = O , (A4)
a=l ~a ~

and all the forces are independent. This can be seen by
summing equation (A2), which results in (A3).

The purpose of this appendix is to show that if the
phenomenological coefficients are related by (A3), then

(A1) requires
2 Q = 0, a=0,...fl . (A5)

Thus only n2 of the original (n+1)2 coefficients are inde-
pendent. de Groot (1952) first noticed the possibility of
(A5) for a simple case, but failed to provide any generaliza-
tion of (A5) or its proof.

Two lemmas will be used in establishing (A5):

Lemma I. The quadratic form

a11 a12 Y1
9(3) = ¥TA¥ = Y1 Y2 <A6>
a12 a22 Y2

127
is positive (q 2 0) for any Y if and only if both

all 2 0 (A7)

and

2

‘ a12

all a22 2 0 . (A8)

Proof. Equation (A6) may be rearranged to give

 

 

a12 2 312 2
q = a ‘Y + ——— Y + a - ——— Y . (A9)
11 1 all 2 22 all 2

 

Thus if for any Y1 and Y2, all 2 0 and alla22 - a12 2 0 then

q 2 0. On the other hand if q 2 0 for any Y, and Y2 = 0:

a
12 2
Y1 + a 11322 " 5‘12 3 0'

then all 2 0, while if Y2) = 0 then a

 

Definition. Let

 

V
36 5 2 add 6:0,..." (A10)
d=l
and
\J
05 = 2 06a 6:0,...n . (All)
d=1

Lemma II.

0:0,...0 (A12)

Proof. For each 0, 0:0,...0, choose

128

§0 = §d = 0 ’ a > v ’
§a = §Y , a=1,...v, y+6, d+5 . (A13)

Then by use of (AlO-ll)

l .
,.__1~ “66 2 (36+C6)‘968 X6
0 = xtS x
.1__;l 1 v . (A14)
7(36+C5)-(266 SE135+Q55-BG-C xY
Use of Lemma I and (A1) requires
V 1 2
066 £138 2 4 (135 + Ca) , 6=0,...n . (A12)
3 =

The final theorem which gives (A5) is now:
Theorem. Given a scalar, invariant, positive,

bilinear function 0, defined by

(Al)

Q
Ill
IIM=I
IL:
w
o
Q

where the X0 are a linearly independent set of vectors
which determine the vectors Ja uniquely according to the

linear relations

TT
~a = 8:0 QanB , a=0,...n , (A2)
then if
\)
) Q08 = 0 , 8:0,...n , (A3)

(A1) requires

Proof.

and therefore

or

By hypothesis

129

(1:0,...71' o

6=0,...n ,

(1:0,00011' 0

(A5)

(A13)

(A14)

(A15)

(A5)

APPENDIX B

PARAMETER ESTIMATION PROGRAMS

l. Subroutine MINIMIZE

The subroutine described here was written in

 

FORTRAN for use with a CDC3600 digital computer. It is
designed to minimize a function of up to ten variables
by choosing conjugate search directions. This assures
that a quadratic function of n variables will be mini-
mized in at most n steps. (If the number of variables
exceeds ten, the program must be redimensioned.) For a
theoretical description see the article by Powell (1964),
"An efficient method for finding the minimum of a func-
tion of several variables without calculating derivatives."
There are three considerations for the use of the
subroutine which must be tailored to the individual pur—
pose.

1. Calling statement. The subroutine may be called

 

from a FORTRAN program by means of the following
statement:
CALL MINIMIZE(X,N,EPS,ENDNORM,ITMAX,IPRINT,SUCCESS,FNORM)

2. Parameters.

 

a. X = A linear array dimensioned for the number

of variables. The program should be called

130

131

with a set of initial guesses for the variables
stored in X. The solution will be returned in X.
N = The number of variables (less than ten).

EPS = A convergence criterion parameter. The
change in each variable from the-last step is
compared with EPS times the current value, and
convergence is assumed if the change is smaller.
ENDNORM = A convergence criterion parameter.

The function value at the current point must

be less than ENDNORM to obtain convergence.
ITMAX = The maximum desired number of iterations.
IPRINT = An option. If IPRINT equals unity, the
program will cause the results to be printed.

If IPRINT is zero, no results will be printed.
SUCCESS = A logical variable to indicate con-
vergence. If success is unity, the process has
converged. If SUCCESS is zero, the method has
failed to converge, and a statement will be
written to indicate the reason for termination.
(SUCCESS must be declared TYPE LOGICAL in the
calling program.)

FNORM = The name of the function subroutine
which calculates the function to be minimized.
(See the next section. FNORM must be declared

EXTERNAL in the calling program.)

132

3. Required subroutines.

 

a. QUADMIN. This is a routine required by MINIMIZE
and a listing is given together with MINIMIZE.

b. FNORM. This is a subprogram in which the function
to be minimized is calculated. It must have the
following form:

FUNCTION FNORM(X,N)

DIMENSION X (N)

(any necessary calculations)

FNORM =- f(X(1), X(2), ..., X(N))
RETURN
END

where f is the function to be minimized, and

X and N have their previously indicated meanings.

 

Listing of MINIMIZE and QUADMIN

 

134

SUBROUTINE MINIMIZE(XQNOEPSIQEPSZQITMAXQIPRINTQSUCCESSQFNORM)
DIMENSION X(N)9XO(10)9Y(10)9R(10910)
COMMON/MIN/LASTNORMoKOUNT

EXTERNAL FNORM

TYRE REAL NORM. LASTNORM

TYRE LOGICAL SUCCESS

IF(N.GT¢10)GO TO 5000

ITER=O

KOUNT=O

D01 I=19N

X0(I)=X(I)

P(IQI)=OoI*XO(I)
IE(X0(I)¢LT¢(1¢OE-7))P(I:I)=OoOI

L=I+l
DO 1 J = LoN
1 P(IQJ)=P(J91) =OoO

LASTNORM = FNORM(X9N)
KOUNT=KOUNT+1
NM=N~1
IF(IRRINT)PRINT lOOoLASTNORMqX
100 FORMATIIHIQ*THE INITIAL VALUES ARE*O//O5XQ*NOQM*QIOX9*XO(1)OOQOOXO(N)*9
I(N)*o//09E1506o/1(15X08E1506))
IF(IRRINT)RRINT 110
110 FORMAT(1H60*ITER INC*95X9*NORM*9lOXo*X(1)oooooX(N)*o//)
1000 ITER =ITER+1
IF(ITER.GT¢ITMAX)GO TO 3000
95 DELTA=1.0E-IOO
M=O
F1=LASTNORM
DO 2000 I=IQN
DO 2 J=IoN
2 Y(J)=P(J9I)
CALL QUADMIN(X CY ONORMQN gFNORM)
IF(IRRINT)PRINT lOloITERoIqNORMcx
1(31 FORMAT(2I598E15.60/o(25X97E1566))
IF((LASTNORM-NORM).GE.DELTA)3,4

3 M=I
DELTA 8 LASTNORM-NORM
4 LASTNORM = NORM

2000 CONTINUE
2004 F2=NORM
IF(ITER.GT.N)15916

155 IF(NORM.GT¢EPSZ)16¢I7
1‘7 IF(ABS(FI-F2)oGToABS(ERSl*F2))16919
1‘? DO 18 I=loN
IF(ABS(X(I)-XO(I))oGToABS(ERSI*X(I)))16918
153 CONTINUE
GO TO 4000
15> DO 5 =10N

£5 Y(I)=2.0*X(I)-XO(I)

135

F3=FNORM(Y9N)

KOUNT=KOUNT+I

IF((F3.GEoFl)o0Ro(((F1-2.0*F2+F3)*(F1-F2-DELTA)**
12).GE.(DELTA*((F1-F2)**2)/2.0)))6o7

6 DO 8 I=1¢N
8 XO(I)=X(I)
GO TO 1000
7 DO 9 I=19N
9 Y(I) =X(I)-XO(I)

CALL QUADMIN(X0YoNORMoNoFNORM)
DO 10 I=IQN
IO XO(I)=X(I)
DO 11 I=M9NM
DO 11 J=10N
11 P(JoI)=P(JQI+1)
DO 12 I=19N
12 P(ION)=Y(I)
LASTNORM = NORM
GO TO 1000
3000 PRINT 102
102 FORMAT(1H49*THE MAXIMUM NUMBER OF ITERATIONS HAS BEEN EXCEEDED*)
SUCCESS =0
PRINT 5004gKOUNT
RETURN
4000 PRINT 1039ITER
103 FORMAT(1H49*THE PROCESS HAS CONVERGED IN*91693X9*ITERATIONS*)
SUCCESS =1
PRINT 50049KOUNT
5004 FORMAT‘lH-O*THE NUMBER OF FUNCTIONAL EVALUATIONS WAS*,110)
RETURN
5000 PRINT 5001
55001 FORMATC1H40*MORE THAN 10 VARIABLES. PLEASE REDIMENSION.*)
SUCCE55=0
RETURN
END
SUBROUTINE QUADMIN(X9P.NORM9N9FNORMI
DIMENSION PHII3)9VT(3)0X(N)9P(N)
COMMON/MIN/LASTNORMQKOUNT
TYPE INTEGER UPPER
TYPE REAL NORM vLASTTo LASTNORM
DO 9 I=loN

9 X(I) =X(I) +P(I)
LASTT = 1.0
T3000
ITER =0

1(3 ITER = ITER + l
NORM = FNORM(X¢N)

KOUNT=KOUNT+1
IF((ABS(T-LASTT)oGTo(001*ABS(T))oANDoITERoLEoEO)oORo<ITERoEQoB>I
111912

11
13

14

15

16
17

18
19
20

21
22

23
24

136

IFIITERoEOoIII3oI4
VT(11=OOO

VTI3) I“1.0

PHIII) =LASTNORM

PHII3) =NORM
IF(PHI(I)oGToPHI(3))10?
T=200
LOWER=I

MID=3

UPPER=2

K=2

GO TO 1000

T=-l.0
LOWER=2

MID=1

UPPER=3

K=2

GO TO 1000

PHI‘II)=NORM

XW =VT‘21‘VTI3)
XX=VT(3)‘VT(I)

XY =VT(1)-VT(2)

XW =-(PHI(1)*XW+PHI(2)*XX+PHI(3)*XY)/(XW*XX*XY)
xx=(PHI(11-pHI(2))/XY-XW*(VTII)+VT(2)I
LASTT =T

IFIXWOGTQOQO)ISQI6
T=-XX/(2.0*XW1

GO TO 19
IF(PHI(UPPER)¢GT.PHI(LOWER))17918
T:3.O*VT(LOWER)‘2.O*VT(MID)
GO TO 19
T=3.0*VT(UPPER)‘2.O*VT(MID)
IF(T.GT.VT(UPPER1)20921
I=LOWER
LOWER =MID

MID = UPPER

UPPER 31

K=UPPER

GO TO 1000
IFCT.LT.VT(LOWER))22923
I=UPPER
UPPER =MID

MID =LOWER
LOWER=I

K=LOWER

GO TO 1000
IE‘T.GT¢VT(MID))24025
I=LOWER

LOWER =MID

MID=I

25

1000

1001

12

137
K=MID
GO TO 1000
I =UPPER
UPPER =MID
MID=I
K=MID
II=K
VTIK)=T
DO 1001 J=19N
X(J)=X(J)+(T-LASTT)*P(J)
GO TO 10
IF(NORM.LE.LASTNORM)RETURN
NORM=LASTNORM
DO 7 I=10N
X(I)=X(I)‘LASTT*P(I)
RETURN
END

138

2. Program FOTOFIT

 

The fringe shape analysis described in section 4
of Chapter IV is performed by the program FOTOFIT and its
associated subroutine FRNGFIT and FRNGSTAT. FOTOFIT is
an executive routine which handles the reading of data,
transformation of data, and calling of auxiliary programs.
The actual calculations of fringe shape and weight func-
tions are done in FRNGFIT while the estimates of standard
deviation are made from the moment matrix calculated by
the subroutine FRNGSTAT. Of course the minimization of
the difference between the calculated and experimental
fringe shapes is done by the routine MINIMIZE which must
accompany these routines.

The data input format is given by FORMAT statement
3 in the listing that follows. The first card of the data
deck contains an identification code, the number of data
cards to be read, the number of derivatives to be estimated,

and the parameters A, b1, K and u1 as specified in equa-

TI
tion (4.6) of Chapter IV. The remaining data cards are
read according to FORMAT statement 2 and each card contains
one (xi, zi) point. As many fringes as desired may be

analyzed at one time. The only requirement is that a blank

card be placed at the end of the data deck.

Listing of FO'I‘OFIT, FRNGFIT, and FRNGSTAT

1000

10

999

140
PROGRAM FOTOFIT
DIMENSION X(lO)
COMMON/TRNS/Z(29100)oNOPT¢BloA9TK9UT1oWTIlOO)9F(100)9TEMP9FOTOLBL
TYPE LOGICAL SUCCESS
EXTERNAL FRNGFITQWTFIT
READ 3OFOTOLBLONOPTONOAQBIOTKQUTI
F0RMAT(A802I594E15.6)
IF(NOPT¢E000)GO TO 999
D0 1 I=IQNOPT
READ 292(19I19Z(29I)
FORMAT(F6¢39F5.3)
2(19I13200*(Z(10I)+0035)/404
CONTINUE
DO 10 I'IQN
X(I13000
CALL MINIMIZE(X0N00000101.00259IQSUCCESSQFRNGFIT)
CALL MINIMIZEIXQNQOOOOOI01000501IOSUCCESSQWTFIT)
CALL FRNGSTAT(X9N)
GO TO 1000
CONTINUE
END
FUNCTION FRNGFIT(XQN)
DIMENSION X(N)
COMMON/TRNS/Z(29100)oNOPToBlvoTKoUTIoWT(IOO)0F(100)9TEMP9FOTOLBL
DO 1 1819NOPT
WTIII=100
LUSE=1
GO TO 2
ENTRY WTFIT
LUSE=0
CONTINUE
M=N-1
DO 3 I'IQNOPT
YP=Z(IQI)+BI
YL=YP‘20*B1
TP=TK*(YP**2-1001/40+YP/20
TM=TK*(YL**2-1001/40+YL/20
IFILUSE)GO TO 5
TEMP=A*UT1*TK+BI
DO 4 J31 QM
TEMPITEMP+A*XIJI*(J+I,*(Tp**J*(TK*YP+IOI“TM**J*(TK*YL+Io)1/20
WTIII=IoO/(100+TEMP**2)
CONTINUE
DN=UT1*(TP*TM)
DO 6 K'IOM
DN=DN+X(K)*(TP**(K+1)“TM**(K+1)1
FCI)=X(N)-Z(29I)+A*DN
TEMP=000
DO 7 I=19NOPT
TEMP=TEMP+WT(I)*F(I)**2

3

2102

3000

2106
2101

3001

3002

301
3003

211

141
FRNGFIT=TEMP

RETURN

END

SUBROUTINE FRNGSTAT(X9N)
DIMENSION X(N)¢D(10010)oS(10v100)

COMMON/TRNS/ZCZQ1001CNOPTQBI.AoTKvUTIoWT(IOO)OFIIOOIOTEMPQFOTOLBL

M=N-1

DO 1 I319NOPT

S(NQI)=Io

DO 1 J=19M

YP=Z(IOI)+BI

YL=YP-20*BI
TP=TK*(YP**2-1oO)/4o+YP/2.
TM=TK*(YL**2‘IoO’/4o+YL/20
S(JoI13(TP**(J+11-TM**(J+1))*A
DO 2 K=IoN

DO 2 J'ION

D(K9J)‘000

DO 2 ISIQNOPT
DIKOJ)8S(KQI,*WTII)*S(JOI)+D(KOJ,
CALL INVERSE(D¢N9N010E~7QDET010010)
IF(DET0E00000)PRINT 3
FORMAT11H6¢*MATRIX IS SINGULAR*)
PRINT 21020(IOPQIOP=10N)
FORMAT(/25X0*SIGMA MATRIX*/169121101
PHI=NOPT~N

BV=TEMP/PHI

D0 3000 I=1QN

DO 3000 J=IoN

D‘IQJ)=D(IOJ,*BV

DO 2106 J=10N

PRINT 21019J0(D(JOI)OI=19N)
FORMAT(* *9I2913E1002)

DO 3001 I’IQN

S‘IOII'SORT(D(IOI))

PRINT 30020FOTOLBL

FORMAT(1H~¢*FINAL RESULTS FOR RUN*9A8q//9*PARAMETER NUMBER*95X0*ES

ITIMATE*913X0*STDo DEV.*)

DO 301 I819N

PRINT 3003910X(I)9$(191’

FORMAT(IH 95X.I5c5X9E150605X0E1503)
BV=SQRT(BV)

PRINT ZIIIBVOPHI

FORMAT(* STANDARD ERROR OF Y *9E15o69*
1F8.2)

PRINT 90769(19F(I)91=19NOPT)

9076 FORMAT(4(4H FO(13q3H) =E15.7.5X))

RETURN
END
C6H12294 29 5 2.394998E+5 2.272727E-1

-1.4001o504

DEGREES OF FREEDOM *9

1.200610E-3

*1.8500C

-1.3001o257
-1o2001o062
-1o10000894
~1.0000.74S
-0.9000o619
-0.8000o514
-O.7000o405
—0.6000o333
-Oo3000.104
-0.20000034
-0.1000o014
00.0000o000
0010000009
0.20000025
0.3000o074
0.4000o159
0.5000o225
0.6000o324
0.7000o460
0.8000o569
009000.691
1o0000o845
1.10010046
1.2001o215
1o3001o524
01.4001o834

142

143

3. Program EXPCURVE

 

EXPCURVE is a program similar in nature to the
FOTOFIT program and used for the time constant analysis
of pure thermal diffusion experiments. Again the program
is simply an executive routine which reads data and calls
on the auxiliary routines EXPO and EXPOSTAT. The expo-
nential functions are calculated together with weighting
functions in EXPO while the subroutine EXPOSTAT handles
the statistical calculations.

The data input is as follows: The first card to
be read contains the number of data points in the time
series to be analyzed. The second card contains the ini-
tial estimates of the parameters a, b, and c in equation
(5.2) of Chapter IV. Each subsequent card is punched with
one (di' ti) data point which is to be read according to
FORMAT statement l03. As listed the routines analyze only
one set of data at a time, but they may be easily repro-

grammed to handle multiple data decks.

Listing of EXPCURVE, EXPO, and EXPOSTAT

101

102

103

145

PROGRAM EXPCURVE

DIMENSION X13)
COMMON/TRNS/ZI2QIOO)ONOPTOTEMPOWT(100)
EXTERNAL EXPOQWTEXPO

TYPE LOGICAL SUCCESS

READ IOIQNOPT

FORMAT(IZ)

READ 1029(X(I)91=103)
FORMATIF6029E4039F502)

DO 1 I=IQNOPT

READ 103921201)!Z(101)
FORMAT(F7029F502)

CALL MINIMIZEIXOBOOOOOOI0005910.IQSUCCESSOEXPO)
CALL MINIMIZE1X930000001900591091OSUCCESSQWTEXPO)
CALL EXPOSTAT(X9N)

END

FUNCTION EXPO(X9N)

DIMENSION X(N)$F(1001
COMMON/TRNS/Z129100)QNOPToTEMPQWTI100)
DO I I=10NOPT

WTIII=100

LUSE=1

GO TO 2

ENTRY WTEXPO

LUSE=O

CONTINUE

DO 3 I‘IQNOPT

S=X(1)*EXP(-X12)*Z(EOI))

IFILUSE1GO TO 4
WTII1310/(Io+(X(2)*S)**2)
F11)=Z(1911-$-X(3)

CONTINUE

TEMP=OOO

DO 5 I=IoNOPT

TEMP=TEMP+WT(I)*E(I)**2

EXPO=TEMP

RETURN

END

SUBROUTINE EXPOSTAT(X9N)

DIMENSION D(30319X(N)OS(391001
COMMON/TRNS/Z129100)QNOPTcTEMPoWT(100)
DO 1 I=IQNOPT
51101)=-EXP(-X(2)*Z(ZQI))
5(291)=Z(2.I)*X(1)*EXP(-X(2)*Z(2.I1)
$139113-100

DO 2 K=103

DO 2 J=193

D(K9J)=Oo0

DO 2 1:1.NOPT
D(K9J)=S(K01)*WT(I)*S(J91)+D(K¢J)

146

CALL INVERSE(D.3.3.1oE-7.DET.3o3)
IF(DET.E0.0.0)PRINT 3
3 FORMAT(1H6.*MATRIX IS SINGULAR*)
PRINT 2102.(IOP.IOP=1¢3)
2102 FORMAT(/25Xo*SIGMA MATR1X*/16o12110)
PHI=NOPT-3
BV=TEMP/RHI
DO 3000 1:1.3
DO 3000 J=103
3000 D(IoJ)=D(IoJ)*BV
DO 2106 J=1.3
2106 PRINT 21019J0(D(J9I10I=193)
2101 FORMAT<* *012913E1002)
DO 3001 1:103
3001 5(1.1)=SORT(D(IcI))
PRINT 3002
3002 FORMAT(1H-.*FINAL RESULTS*9/.*PARAMETER NUMBER*o5X.*ESTIMATE*.5X.*
ISTD. DEV.*)
DO 301 1:193
301 PRINT 3003oIvXII)oS(IQI)
3003 FORMAT(1H .SXo15.5x9E15.6.5X.E15.3)
BV=SORT(BV)
PRINT ZIIOBVQPHI
211 FORMAT(* STANDARD ERROR 0F Y *.E15.6o* DEGREES OF FREEDOM *.
IFS-2)
RETURN
END
33 F6037.06833.032-0.2501.7500.9766934.950
-OS.50.01610.07
0025.0005.71
0030.0005.97
0035.0006.22
0040.0006.54
0045.0006.75
0050.0007.13
0055.0007.30
0060.0007.63
0065.0007.90
0071.0008.16
0075.0008.39
0080.0008.45
0090.0008.66
0100.0008.78
0110.0009.11
0120.0009.31
0130.0009.41
0140.0009.46
0150.0009.74
0160.0009.58
0170.0009.82

DISC-0009.72
019000009085
020000010007
0212.0009o99
0222:0009-82
0240-0009o84
0252.0009o97
0266.0009.96
0280.0010o01
030000009099
0314.0010007
0364-0010007

147

APPENDIX C

VELOCITY AND TEMPERATURE PROFILES FOR THE

CONVECTIVE HEAT LOSS PROBLEM

l. Velocitnyrofiles

 

The detailed solutions for the vertical and radial
velocity components in a pure thermal diffusion experiment
with radial heat loss have been given in Chapter III,
section 5. It is shown there that the velocities are

determined by the potential ¢ according to

-1
_ M 1 30
V10x "' [doc + F] P ’5? “3'“
-l
_ _ M 1 34>
VlOz " [doc + 17'] a? a? “-2)

where r and 2 have been scaled to unit length and where,

by equation (III.5.16),

-1
M _ w
[doo + T] " M

The potential ¢ is given in equation (III.5.62) as

2 2 N M
O = r(r-l)z(z-l) r(r-l) + 2 2 C
i=0 j=l

..zl
1]

148

exp [A(1/2 - 2)] (C.3)

sin1rjr (C.4)

149

The coefficients Cij have been determined by the minimiza-
tion of the integral in equation (III.5.65) using the
MINIMIZE routine (see Appendix B). It_is clear that
their values depend only on the magnitude of the quantity
A. The numerical values determined with A equal to 0.01

are
0.12007 - 0.00597 0.00105

= - 0.17734 - 0.00409 0.00197 (C.5)
0.17785 - 0.00396 0.00201

Having determined these coefficients, one can

easily calculate vlOr and v10 for any r and z. The

2
vertical velocity lez has been plotted versus r for 2
equal to 0.5 and is shown in Figure 1. For other values
of z the vertical velocity has the same shape but de-
creases in magnitude almost parabolically from the center
towards the horizontal plates.

The radial velocity vlOr has been plotted versus
r for z = 0.25 and the results are shown in Figure 2.
Again the shape of vlOr remains essentially the same for

all values of z. The magnitude of v r decreases

10
as 2 approaches 0.5, where ler changes sign, and then
increases in magnitude.

It should be recognized that the signs on the
velocities as they have been plotted are not absolute,

but should be multiplied by the sign of the heat loss

coefficient. Thus if there is a heat loss, the fluid

150

102

 

 

 

Figure l.-&Vertical velocity as a function of r at
the center of the cell (2 = 0.5).

vlOr

0.25

0.20

0.15

0.05

151

 

Figure 2.--Radial velocity as a function of r at z = 0.25.

152

near the wall will be relatively more dense than that
in the center and the velocity near the wall will be

downward.

2. Temperature Profiles

 

The temperature correction, T10 (r,z), for the
convective heat loss problem is determined according to

equation (III.5.73) by

00

T10 (r,z) = nil Tn(r) sin (nnz) (C.6)

where

Tn(r) = ClI0 (nnr) + Tnp (r) ’ (C.7)

and the coefficient C1 is given by equation (III.5.83).
The zeroth order modified Bessel function-has been denoted

by I0 (mrr). The particular solution, '1' (r), is given by

np

Tnp (r) = 220 Gng’r2 (C.8)
and the coefficients Gn2 are given by equations (III.5.77-8).

The temperature correction depends on the phenome-
nological coefficients for the system through the three
parameters A, pl, and p2. The coefficient A is again taken
to be 0.01 while pl is taken to be 100 and p2 to be 500.
These are values typical to common liquids. The shape of

T is also affected by the value of the quantity

lO

153

(Tr - Tm + 1/2 AT) which enters through the boundary
conditions. The sign of this quantity determines
whether there is an overall heat loss or heat gain
(- and + respectively).

The temperature correction has been plotted
versus r at the center of the cell (2 = 0.5) under the
three conditions (Tr - Tm + l/2 AT) = - 5.0, 0.0, and
5.0 in the Figures 3, 4, and 5 respectively.

It is interesting to note that the lepe of T10
versus r has the same sign as (Tr - Tm + l/2 AT) as r
approaches l.0 indicating heat loss or gain as the case
may be. The temperature corrections appear to be largest
near the center of the cell (2 = 0.5) and to decrease in

either direction towards the horizontal plates.

10

-6.0

-700

-800

 

154

Figure 3.—-Temperature correction versus r at z = 0.5

10

155

 

 

Figure 4.--Temperature correction versus r at z = 0.5
for (Tr - Tm + 1/2 AT) = 0.0.

10

-405 E

 

Figure 5.--Temperature correction versus r at z

156

for (Tr - Tm + l/2 AT)

5.0.

0.5