^FLtJID FLOW IN A DIFFUSELY
REFLECTING RECTANGULAR CHANNEL
AT HIGH SPEED IN THE FREE
MOLECULE REGIME AT ZERO ANGLE
OF ATTACK
B.S., M.E., University of Colorado, 1983
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Department of Mechanical Engineering
by
Dewey E. an
Master of Science
1985
This thesis for the Master of Science degree by
Dewey E. Widman
has been approved for the
Department of
Mechanical Engineering
by
Ali U. Ansari
iii
Widman, Dewey E. (M.S., Mechanical Engineering)
Fluid Flow in a Diffusely Reflecting Rectangular Channel at High
Speed in the Free Molecule Regime at Zero Angle of Attack
Thesis directed by Assistant Professor Ali Ansari
Approximate closed form solutions are found for the free
molecule integral equations that represent the transfer of mass,
energy, and momentum in a diffusely reflecting rectangular channel
at zero angle of attack at high speed. A separable kernel of the
exponential form is used to approximate the actual kernel. Based on
this, closed form solutions are obtained in a form which makes para
metric dependencies very clear, in addition to allowing one to do
handbook type calculations. Approach to the external flow (single
plate) limit is studied and conclusions are drawn with respect to
this limit as a function of various parameters. Plots of heat
transfer and lift are shown for a large range of plate temperatures
and channel dimensions
iv
ACKNOWLEDGMENTS
I would like to give my sincere thanks to Dr. A. Ansari, who
provided continual help and support in the formulation and writing
of this thesis, and to my other committee members, Dr. K. Ortega and
Dr. J. Trapp, whose comments and suggestions gave me a better
insight into the nature of the problem.
CONTENTS
ABSTRACT......................................................... iii
ACKNOWLEDGMENTS.................................................. iv
LIST OF FIGURES.................................................. vii
NOMENCLATURE.................................................... viii
CHAPTER
I. INTRODUCTION................................................ 1
II. STUDY OF A RECTANGULAR CHANNEL.............................. 6
Number Flux............................................. 6
Energy Equations.......................................... 8
Conservation of Momentum................................. 11
Calculation of Drag.................................... 12
Calculation of Lift.................................... 12
III. AN APPROXIMATE SOLUTION IN THE HIGHSPEED LIMIT............ 14
Conservation of Mass..................................... 14
Conservation of Energy................................... 15
Conservation of Momentum................................. 19
IV. RESULTS AND CONCLUSIONS.................................... 21
Number Flux.............................................. 21
Energy and Heat Transfer................................. 21
Momentum................................................. 23
Calculation of Lift.................................... 23
Conclusions............................................ 27
vi
BIBLIOGRAPHY..................................................... 29
APPENDIX
1. DERIVATION OF C AND c2..................................... 30
2. THE DERIVATION AND SOLUTION TO EQUATIONS (3.5)........... 32
3. DERIVATION OF THE CONSTANTS IN EQUATION (3.7)............ 35
4. COMPUTER PROGRAMS.......................................... 40
vii
FIGURES
Figure
1 1 Rectangular Channel.......................................... 3
3.1 ^ext at sPeed......................................... 16
3.2 E at High Speed............................................. 18
4.1 Number Flux at High Speed.................................. 22
4.2(a) Q for T = .1.............................................. 24
1 1
4.2(b) Q for T = 2................................................. 25
1 1
4.2(c) Q for T = 5............................................ 26
1 1
4.3 Lift on Plate 1............................................. 28
NOMENCLATURE
Area
Constant used in approximating the probability kernel
Specific heat at constant pressure
Specific heat at constant volume
Energy flux (energy/area time)
Separation distance of flat plates
Number flux intensity (#/area time solid angle)
Momentum flux intensity (momentum/area time solid angle)
Energy flux intensity (energy/area time solid angle)
Heighttolength ratio for the channel H/L
16 o
Boltzmann constant (1.38 10 ergs/ k)
Probability kernel
Knudsen number
Length of the channel
Mass of the molecule
Number flux (#/area time)
Number density
Momentum flux (momentum/area time)
Heat flux (heat flow/area time)
Gas constant
Position vector
Speed ratio (V /V )
o mp
ix
T
Vo
Vmp
x,y,z
a
gf 6
Y
a
X
T
?
1
2
i
r
w
0
1
r
T
n
Temperature
Macroscopic velocity
Most probable thermal speed of a molecule
Cartesian coordinates
Energy accommodation coefficient
Angles
Ratio of specific heats (C /C )
p v
Momentum accommodation coefficient
Mean free path
Tangential stress
(1a)
Subscripts
Upper plate
Lower plate
Incident
Reflected
Wall
Opening
Left
Right
Tangential component
Normal component
Superscripts
Nondimensionalized quantity
CHAPTER I
INTRODUCTION
In studying fluid flows, three regimes may be analyzed:
continuum, transition, and free molecule. These three regimes can
be characterized by different values of a parameter called the Knud
sen number kn = A/L where A is the mean free path of the particles,
and L is a characteristic dimension of the channel. Fluid flows in
which kn<<1 are classified as continuum, kn>>1 as free molecule, and
kn=!l as transition. This paper deals with the study of free mole
cule flow only. In studying free molecule flow, the assumption is
made that the transfer of energy is due only to the collision of
particles against the surface and that particletoparticle colli
sions are rare.
An important application of the problem of free molecule
flow in a rectangular channel is mainly in the area of space vehi
cles and high altitude flight at high or moderately high speeds.
Much effort has been spent in the calculation of drag and heat
transfer for various geometries, but all of these calculations
involve various numerical methods, most notably the Monte Carlo
method, used by authors too numerous to mention, matrix probability
theory (9), and a variational method used by DeMarcus (5). A closed
form solution has only been found for the cylinder and the sphere
problem (3). The main motivation of this paper is to find
2
approximate closed form solutions to the free molecular transfer
equations for a rectangular channel. Also studied is the external
flow limit in which the distance between the plates is allowed to
increase until the flow over each plate is unaffected by the other
plate. This limiting heighttolength ratio, He is found to be a
function of plate temperature and flow speeds. Many types of space
vehicles have configurations which can be modeled either as a rect
angular channel or as a flat plate.
In the solution to the heat transfer between_j?lates, a
closed form solution is found in the form q = AT^ + BT^ + CS^ + D
where the coefficients depend only on H, the ratio of height to
length of the channel and are independent of the temperature of the
plates or the speed (S). This not only simplifies the presentation
of the data considerably and allows one to do handbook work in cal
culating the heat transfer but it also clarifies parametric depen
dence of all the important variables.
In assuming that the z dimension (see fig. 1.1) is infinite,
we assume that there will be no contribution to the transfer of
energy in that direction. By using this type of geometry, we can
study the limit as the separation between the plates becomes large
enough to yield the solution to the problem of external flow. By
external flow, we mean the problem of a single plate. To study this
limit, we define the limiting value of H, viz. at which the
heat transfer for each plate is a specified approximation of the
external flow or single plate heat transfer. This is a func
tion of the plate temperatures (nondimensionalized with respect to
3
00
LENGTH = L
Figure 1.1 Rectangular Channel
4
the ambient temperatures), the speed of the vehicle, and the energy
accommodation coefficient, which measures the ratio of incident
energy to the reflected energy, and will be described in more detail
later.
As previously mentioned, free molecule flow is characterized
by kn>>1. Transfer of energy takes place by particletosurface
collisions, and particletoparticle collisions are assumed to be
rare. The interactions between the particles and the surface will
depend upon the roughness and cleanliness of the surface. We assume
that no particles are absorbed at the plates. To study energy and
conduction, the usual practice is to assume that a certain portion
of the energy and momentum of incident particles is absorbed and the
rest reflected. The parameters that relate this portion of parti
cles or energy are called accommodation coefficients. In this paper
these parameters will be assumed to be constant over the surface in
question. Three different coefficients will be used. An energy
accommodation coefficient, a
a =
e. e
1 r
e ,e
1 w
(1.1)
where e^is the incident energy and e^ is the reflected energy, and
e is the energy the reflected beam would carry if it came into
w
thermal equilibrium with the wall. If all of the energy is
absorbed, the a = 1. We will also define a tangential momentum
accommodation coefficient
a
T
(1.2)
T
where P represents momentum and x refers to the tangential compo
nent. Again i and r refer to incident and reflected.
5
Similarly, we also have a normal momentum accommodation
coefficient
a
n
P.
1
P
r
n n
P.
1
P
w
n n
(1.3)
where P is the normal component of momentum if the particle
n
distribution achieves equilibrium with the wall.
CHAPTER II
STUDY OF A RECTANGULAR CHANNEL
The sketch of a rectangular channel is shown in fig. 1.1.
The flow is in the x direction and the z direction is much greater
than the y or x direction. Before continuing the analysis, all of
the flux terms will be nondimensionalized with the corresponding
ambient equilibrium gas fluxes. These expressions are
mp
n = N /
00 00 2 T7
00 v+1
e = mRN 1
oo 00 4 v it Y1
3mN V2
P =
00
mp
Also, we will define a nondimensional speed S = v /v where v
o mp o
equals the speed of the reference frame; i.e., the flight speed of
the space vehicle. Plate temperatures will be nondimensionalized
with respect to the ambient temperature, T = T/Too. Throughout the
remainder of the text, it will be assumed that the number, energy,
and particle flux terms have all been nondimensionalized with
respect to the above quantities.
Number Flux
By considering the geometry of fig. 1.1 with subscripts 1
and 2 representing the upper and lower walls, we can write the
number flux equations as
n (r ) = n (r ) + /, n (r )k(r r )dA
1 w 1 ^ w. 1 Aw_ 2 w w w w
1 ext 12 212
n (r ) = n (r ) + / n (r )k(r ,r )dA
2 w_ 2 . w_ J A 1 w. w. w. w.
2 ext 2 w. 1 12 1
7
where ngxt is the function that describes the flux of incoming
(ambient) particles from the left and right openings. In writing
this equation we assume that the particles reflect in a diffuse
manner and that there is no accumulation of particles on the sur
face. k(r ,r ,) is a probability factor which relates the proba
bility that a particle leaving the point w' will arrive at w. The
evaluation of n^^ is arrived at by integrating the number flux
intensity coming from the outside i^Cr^,) over the opening
n = / l (r ,0(r ,r ))k(r ,r )dr .
ext 'Aoo ow ow o
o
where the number flux intensity in(r ,0(r ,r )) is found by assuming
o o o w
a Maxwellian distribution of particles,
___2
(vvo)
____ if 2RT
f (r ,v) = N () e
o o o 2irRT
o
here N T are the ambient conditions for number density and
o o
temperature, and v is the bulk velocity. This function is then
o
integrated over all angles and all particle velocities.
. n
i = ff (r v)dr ,dr
o J o o o
The result is
N v 2
_ _ o mp 2 S i *
in(r ,6) = f e"S (1+Sp + /irS e 6( + S ) (1 + erf(S ).
OO _ j/ Â£ o Q ^ U o
2tt
x 2
X
where SQ = Scos(0), erf(x) = 2 Je dx, V is the most probable
m o mpo
particle speed in a stationary Maxwellian gas at temperature T and
V
S is a speed ratio, defined by S = . For details of the
mp
integration see (1)
8
In considering the space vehicle problem at zero angle of
attack, we assume that the flow over the channel will be uniform and
hence n. = n Furthermore, n = n_ = f(x) = f(r )+f_(r ),
12 1 ^ 2 ^ 1 w 2 w
ext ext
where f^ and f^ represent the number flux coming from the left and
right, respectively. Under these assumptions the number flux
problem reduces to a single equation
n(r ) = n (r ) + f n(r ,)k(r ,r )dA
w ext w JA w w' w w
w'
(2.1 )
The function f(x) or (x), which represents the number
flux coming from both the right and left sides can be shown to equal
(1)
2 2
f(x) = f^(x) + f2(x) = 1 ^ (ue s u (1 + erf(su))
2 2
_q I 1 _ i )
+ u'e M (1 erf(Su1)))
u =
(x2) T
2 2
x + Hj
u' = u(1x).
With lengths nondimensional, the ranges for x and y are 0 < x < 1,
and 0 < y < H and < z < + .
At this point the functions n and k(r ,r ,) are still
r ext w w'
dependent upon z. After integrating this dependence out We get for
the final result (1)
2 1
, . , H r n(x') , ,
n(x) = f(x) + J dx'
Z 0
22 2
((x'x)z+h zr
where f(x) = and x' denotes the plate opposite from x.
(2.2)
Energy Equations
Using subscripts 1 and 2 to indicate upper and lower plate,
respectively, and assuming constant energy accommodation coefficient
9
we have,
e. (r ) = e. (r ) + a/, e (r )k(r ,r )dA
, i. w. 1 w. 1A w w' w.
1 1 ext 1 
w2 "2 "2 "2 "1 W2
+ (1ct)/_ e. (r )k(r ,r )
J A w w w.
dA
w2 2 2 2 1
w_
e. (r ) = e. (r ) + af e.(r )k(r ,r )dA +
i_ w t w_ J A 1 w. w. w w.
2 2 ext 2 w. 1121
(1a)/,. e. (r )k(r ,r )dA
J A l w. w_ w w.
w. n 1 12 1
1
where e.(r ), and e (r ) are the incident and reflected energies,
l w r w
6ext^rw^ t^1Â£ externa^ forcing function for the incident energy.
,e
e. = [ i (r e(r ,r ))k(r ,r )dr
l, J A o o ow ow
ext o
0 n
where i is evaluated in a manner similar to i The result is
o o
mRN V T
o mp o 2
. e, _, o S
W6 .3/2 '
2tt~
j so r
is 6 i
4 6 [
15 + 20Sq + 4Sg
f !*e J]
3 [ + erÂ£1Se]]
+ 4~3T
mRT in(r ,0)
2(y1) o o o
where = Scos0, y is the ratio of specific heats, y = c /c^, and
in(r ,0) is given above. In order to relate the reflected energy
o o
flux to the incident, we must use equation (1.1)
e = (1a)e + ae
r i w
where e (r ) is by definition the energy that would be carried by
w w
the equilibrium Maxwellian flux temperature T, and from kinetic
10
theory
w w
,4^11 n(7.............
2 ( y1 ) WWW
e (r ) = mR ) 1' n(r ) T (r )
By using the fact that e (r ) = T (r )n(r ) in nondimensionalized
w w w w w
form and equation (2.2) we obtain
[a e (r )k(r ,r )dA = T fn(r )k(r ,r )dA
J w w w w w, w 2 J w w w_ w
222 212 2 212
= T (n(r )f(r ))
2 wi w
1 1
and likewise for wall 1. After integrating the function e. over
ext
< z < + and 0 < y < 1, we obtain (1)
e. (x) = f(x) + 2 \
i (y+1)
ext
s2 s2 / :
 e  u
fs2 s_
L 2 / TT
S2, 2 ,2.
e (u u' )
(2.3)
2 ,2
3eSu (1+erf(Su) + u,3eSu (1erf(Su
)]
We obtain two simultaneous equations for the incident
energy fluxes.
e. = e. (x)+aT (n(x)f(x))+(1
i_ l . 1
2 ext
a)f Je (x')dxC1^1
0 1 Lx'x) +H J
2 (2.4)
e. (x) = e. (x)+aT_(n(x)f(x))+(1
1 2
1 ext
a)f Je (x')dx'C
0 2 Lx'x) +H
The heat transfer to plate 1 and 2 may now be calculated
q (x) = ct(e. (x)e (x)) = a (e. (x) T n(x)) (2.5)
1 , *1 1
q2(x) = a(e (x) T2n(x))
12
1
q. = /q.(x)dx
0
and likewise for
11
Conservation of Momentum
Again, the equations are derived in the same manner with the
same assumptions, diffuse reflection, etc., but this time we will
consider the normal and tangential components separately. Using n
to denote normal and j, tangential
P. (r ) = P, (r ) + f P (7 )cosg(r ,,7 )k(7 ,,r )d7 (2.6)
1 w l w JA r w' p w' w w' w w
n
ext
w
P. (r ) = P. (r ) +  P (r ,)sin0(r ,,r )k(r ,,r )dr ,(2.7)
l w l w JR r w' w w' w w' w w'
T T w
ext
where P. (r ) = f iP(r r0(r ,r ))cosa(r ,r )k(r ,r )dr
l w JAoo ow ow owe
ext
P. (r ) = / iP(r ,0(r ,r ))sina(r ,r )k(r ,r )dr
l w JR o o ow ow owe
T . O
ext
and the momentum flux intensity is given by (1)
_2
_ MRN T 2 I*
 TTTf e"S [<5se + 2se> * 6 < Â§ + es2 + 2s^i *
2ir"
]
(1 + erfSQ)
where S = Scos(0), and P (r ) = mn(r ) J 2irRT (r ) frem kinetic
0 r o 4 w r w
theory, and 0, a, are the angles between the normals to r^ and the
vectors
r r I and r r
o w I w' w
Once we have calculated the particle flux and reflected
temperature, we can now find the momentum in nondimensionalized
/
2
P = nU )(T (7 ))2
r w r w
form
12
Calculation of Drag
Assuming once again that the reflections are diffuse so that
the reflected tangential momentum flux is zero and, therefore, the
shear stress X (r ) = P. (r ). Using the expression from 2.7 we
w w i,x w
get
(r )=P. (r )+/a P (r )sin6k(r ,r )dA
w_ w i,x w w_ r_ w w w w
11 ext 1 2 2 2 2 1 1
x (r ) = P, (r ) + /A P (r )sin6k(r ,r )dA
W2 W2 lText W2 W1 ri W1 W1 W2 W2
Integrating the expressions for P^^ 00 < z < + 00 and, 0 < y < H, we
get
px (x) = TTÂ¥
ext
S 2 .2. 2
e (u u' ) S
2Srise
2 2
, 3 Su . , ,3 Su'
u e (1erf(Su)) + u' e (1erf(Su
E
where u and u' are defined as before.
,]]
Integrating over z one obtains (1)
2 1
w
 ' /T
(r ) P ' (x) ~ / n(x') / T (x1)
1 1 ext 0 2
(x'x)
22 2
2 ((x'x) +H )
dx' (2.8)
VV  J* f l <*
2 2 ext 0 1 ((x'x) +H )
where n(x) and were previously solved.
Calculation of Lift
The net force acting on the plates will be the difference in
the net force produced by the normal forces.
F_ =f F (r)dA / F (r_ )dA
L A n w, A n w w
system W1 1 1 w2 2 2 2
(2.9)
13
where
F (r ) = P (r ) + f P (r )cos(0) k(r ,r )dr
n. w. n. w JA r w u w. w. w
11 1 ^ 1 w 2 2 2 1 2
ext 2
+ P_ (r )
r,n1 w1
F (r ) = P (r ) + f P (r )cos(q) k(r ,r )dr
n2 2 2 *2 \ ri W1 2 1 1
ext 1
+ P (r )
r ,n w
2 2
If the diffusely reflected flux is in equilibrium, then
P = 2 = \ n( x) /t" .
r,n2 3 3
Again substituting in the expressions for the kernels and the
expressions from (2.7) and the integrations in the y and z direc
tions are carried out, we obtain the following result for the
average lift
2h2 1 r=
?T = ~ /n(x)(/ T (
L 3n > r
system 0 2
Y'(x)j E(u')
K(u)
xy(x) d
where y(x) = v / 2 2 x +H
(x) + x
(2.10)
dx.
integrals of the first and second kind, respectively.
CHAPTER III
AN APPROXIMATE SOLUTION IN THE HIGHSPEED LIMIT
It is the main motivation of this paper to discuss the high
speed limit (i.e., S^7) and to find the conclusions and simplifica
tions that result from this limit. It will be shown that for the
highspeed case, approximate closed form solutions can be obtained
provided a separable kernel is used. In all three casesmass,
energy, and momentumthe integral part of the equation involves an
inseparable kernel; therefore, differentiation of the equation does
not eliminate the integral sign. Choosing a function k(x,x') =
 xX 1 I
ce I I and by successive differentiation, one can eliminate the
integral sign..
Conservation of Mass
In equation (2.2) we replace the kernel, k(x',x) by a
" 3 I x x I
separable kernel c e 1 I, where 3 and c are determined from a
least squares approximation. Before performing successive differ
entiation using Leibnitz' rule, the integral equation is written in
n(x) = n ,(x) + [n(x')ce lX X I dx' + fn(x')ce X I dx' (3.1)
ext 1 J
0 x
Differentiating twice, we obtain (8)
n''(x) + 3(2cg)n(x) = n (x)^2]!^ (x) (3.2)
The success of this method depends on how well the kernel
can be approximated by the exponential. Experience has shown that
15
by choosing nine data points and using a leastsquares method to
calculate 3 and c, a good approximation to the kernel can be found.
For the highspeed case (8=7 or higher), nex^.(x) goes to one
(see fig. 3.1), which considerably simplifies the solution to (2.1).
The solution to a differential equation of the form (3.2) for the
number flux is
n(x) = c^cos (/~kx) + c^sinf/ kx) 3 /k (3.3)
where k = 5(2c3) and c^ and c^ are arbitrary constants. By putting
this expression for n(x) back into the integral equation and evalu
ating it at the endpoints, one will obtain two equations which can be
solved for c^ and c^. The resulting expressions are
RHS(X2 + sin (vHk) + )
C1 (X2+sin(V k)X4) (1X4)X2( lX^cost/ k)+X3) (3'4)
c (1X cos(/~kx) + X_)
1 1 3
C2 (X2 + sin^ kx) + X^)
X1 = (^)(e ^(Pcosf/k) + Vrksin(/_k) + 3)
X2 = (r^)(e ^(3sin(i/~~k) / kcosfi/k)) + / k) 2p
X3 = (7) (3cos(v" k) + /""ksinf/ k) 2p 1 TO
X4 = (^) (5 sinfi^lc) kcos(^ k) + / ke ^ )
RHS = t + Bf + Pc(B1) k k
For a derivation of these constants see Appendix 1.
Conservation of Energy
Once again we approximate the kernel by an exponential and
eliminate the integral sign by successive differentiation. However,
this time we obtain two fourthordered differential equations since
16
Figure 3.1 Next at High Speed
17
we started with two coupled integral equations (equations 2.4),
e^V(x) + (2e2)e1(x) + ( g4 (2 Â£cg)2 ^ (x) = F^x) 
2g2F21 (x) + g4F2(x) 2?cg3F1(x) + 34(A2) + 2^^.
e1V(x) + (2g2)e(x) + ( g4(2Â£cg)2)e (x) = F^V(x) 
12 '
2g2F^1 (x) + g^tx) 2CcB3F1(x) + g4(A1) + 2?cg3A2.
where F (x) = off (c cos(/Tx) + c sin (/ kx))
1 / z 1 / z 1 z
. 2(y1)_2 T T ,g2,
and ^1,2= 1 + vT~s ~a 1,2 a 1,2 (r}
? = (1a)
The solutions to these two fourthorder differential
(3.5)
equations are
4 m. x
e(x) = Â£ a.e + X.cos/~kx + Xsin/""kx + X.
1
i=1
1
1
(3.6)
4 m.x
e2(x) = 2 b^e 1 + Y.jCOs(/~"kx) + Y2sin(/ kx) + Y^
i=1
2(y1)
(provided that goes 1 + ^ ^^. See fig. 3.2, which is
normalized with respect to this quantity).
(VTi)ci
"here X1 TlifrifTt *2 ~ okhV k) u = 1'> mi *7 b22Sc8
2gc(1C)(gT + 2?cT )
x = 21_ + _Rg_
3 (g2c)(g2(2?c)2) P_2?C
For a derivation of the fourthordered differential equations and
the solutions to these equations see Appendix 2. The expressions
for b_^ and Y^, Y2# and Y^ have the nice property that b^ = a^, and
b = a b = a b = a., and that by switching T and T in the
Z Z J O i *i I Z
expressions for X^, X2, and X3 one will get Y1, Y2, and Y3<
18
Figure 3.2 Eext at High Speed
Egxt is normalized with respect to 1 +
2(Y1) 2
(Y+1)
19
In order to evaluate the heat transfer q, equations (2.5)
were used, except that in the expressions for incident energy flux,
e^, the dependence of equation (3.6] upon temperature was separated
out so that a final expression for the heat transfer could be writ
ten in terms of coefficients that are independent of temperature
(and depend on H only).
2 = A (H)T + B(H)T2 + C(I)S2 + D(H) (3.7)
In order to separate out this dependence on temperature, one must
solve the four equations in the four unknowns (i.e., a^) by hand. A
derivation of these coefficients is given in Appendix 3. Since
2
there is a dependence on S in the expression for the incident
energy flux e^(r^) (the speed ratio),, this can also be separated
out. Finally, there are terms left over in e.(r ) which do not
depend on T or S, and these terms were grouped as a constant D(H).
Writing the equation in this way allows one to isolate the effect of
each parameter on the heat transfer.
Conservation of Momentum
For the highspeed case (8=7), the P term in equation
Text
(2.8) equals (1)
4S
3i/ if
The average shear stress is the integral of equation(2.8]
X = / x (x) dx
w J w
Because of the fact that the integral term in equation (2.8) is odd
symmetric about the point x = .5
11 / T (x')(x'x)
/ / n(x')  dx'dx = 0
o o ((x'xr+irr
20
Equation (2.8) now becomes trivial and
T T 4S
w = w = i (3.8)
1s1 2 v i
The calculation of the lift is more complicated. No simple
expression is possible for the highspeed case, and equation (2.10)
must be integrated numerically. A plot of lift versus H is given in
fig. 4.3.
CHAPTER IV
RESULTS AND CONCLUSIONS
Number Flux
The values of number flux (equation 3.2) for high speed
using the approximate kernel are shown in fig. 4.1. These values
confirm the results of (1) in which an exact kernel was used and
spectral theory of linear operators^ was used to solve the integral
equation for the number flux problem. As H approaches infinity, the
kernel in equation (2.2) goes to zero and the number flux approaches
one. By studying the number flux distribution at different values
of H, we can determine, within any desired accuracy, how closely we
have approached the external flow limit (i.e., n(x) = 1, see fig.
4.1). It should also be pointed out that for small H a singularity
in the kernel is approached, and, therefore, the numerical solutions
also become less accurate at small H (i.e., H<0.8).
Energy and Heat Transfer
As mentioned in Chapter II, in evaluating the heat transfer
(equations 2.5) the energy was expressed in terms of coefficients
that were independent of temperature and depended only on H. (For
q one would only switch T and T in the equation for q .) Once
^ \ Â£ 1
Vor a more complete discussion of the spectral theory
solution techniques see (6).
22
l
Figure 4.1 Number Flux at High Speed
23
these constants were evaluated for different values of H, q^ could
be calculated for different combinations of and Figure 4.2
contains the results for different values of temperature.
A comparison of the values for and q^ was made between
the spectral theory results (see 1) and the approximate kernel
method, and very close agreement was found.
From figures 4.2 (a)(c), we see that for equal plate
temperatures, the heat transfer decreases as H increases indicating
that reflections from the opposite have less effect. For unequal
plate temperatures, one plate will act as a source and the other a
sink; however, this temperature difference must be large enough to
overcome the incident energy of the incoming flux of ambient parti
cles; otherwise the heat transfer will behave as if the plate tem
peratures were the same. For example, in fig. 4.2 (b), for T^ =2,
and T = .1 through = 2, the temperature difference is not great
enough to overcome the effect of the incoming flux of particles and
q^ behaves as if = T2" For Tj = ^ and T = 1 through = 2,
temperature difference is large enough to make the opposite plate
behave as a sink so that as H increases, q^ increases.
Momentum
Calculation of Lift
It has been pointed out (1) that by considering the physics
of the problem, it can be seen that the average angle between the
trajectory of the incident coming from the opposite plate and the
normal to the plate decreases as H increases, and that, therefore,
the component of momentum increases and F^ is expected to increase.
24
II
Figure 4.2(a)
Q for T = .
1
25
i
Figure 4.2(b) Q for =
2
26
I
Figure 4.2(c) Q for T =
5
27
This is clearly evident by looking at fig. 4.3 (where positive lift
is defined in the direction of the normal on the lower plate
directed toward the upper plate). As H increases, the absolute
value of the lift is increased until an external limit is reached at
H5. Notice that this lift is an unusual phenomenon, since it
occurs at zero angle of attack. However, when comparing this lift
to the tangential shear stress at high speed (equation (3.8)), we
see that it has a very small order of magnitude.
Conclusions
This paper has shown that by approximating the kernel with
an exponential function, the integral equations resulting form the
flow of particles in a diffusely reflecting rectangular channel in
the free molecule regime can be reduced to a differential equation
which can be solved to give approximate closed form solutions for
the mass, momentum, and energy equations. Only the highspeed case
was considered which allowed one to set the external forcing func
tion equal to one and made the resulting differential equation sim
ple enough to be solved by hand. A more general study in which the
external forcing function is left in general terms results in an
exceedingly complex complementary solution to the differential
equation and was not considered here.
The resulting closed form expression for the energy equation
allowed the equation for heat transfer between the plates to be
expressed in terms of coefficients that were independent of tempera
ture. This allows one to do heat transfer calculations with rela
tive ease with a wide variety of plate temperature combinations.
28
T2=.1
I
Figure 4.3 Lift on Plate 1
BIBLIOGRAPHY
1) Ansari, A., "Free Molecular Drag and Heat Transfer in a Dif
fusely Reflecting Rectangular Channel at Zero Angle of Attack,"
Ph.D. Dissertation, University of Florida, 1973.
2) Ansari, A., Irey, R. K., "Free Molecule Flow in Rectangular
Channels at Arbitrary Speeds," International Symposium on Rar
efied Gas Dynamics, Ninth Symposium, DFVLR Press, Goettingen,
Germany, 1974.
3) Chahine, M. T., Rarefied Gas Dynamics, Second Symposium, Supple
ment 1, Academic Press, New York, 1961.
4) Chien, F., Rarefied Gas Dynamics, Sixth Symposium, vol. 1,
Academic Press, New York, 1969.
5) DeMarcus, W. C., Rarefied Gas Dynamics, Second Symposium,
Supplement 1, Academic Press, New York, 1961.
6) Friedman, B., Principles and Techniques of Applied Mathematics,
Wiley, New York, 1956.
7) Klavins, A., Rarefied Gas Dynamics, vol. 51, Part 1. Progress
in Astronautics and Aeronautics, Princeton University, Prince
ton, New Jersey, 1975.
8) Ozisik, M. N., Radiative Transfer and Interactions with Conduc
tion and Convection> Wiley, New York, 1973.
Wu, Y., Rarefied Gas Dynamics, Fourth Symposium, vol. 1,
Academic Press, New York, 1969.
9)
APPENDIX 1
DERIVATION OF c1 AND C
By subsituting equation (3.3)
n(x) = c^cos ^x + c^sin / k x
into the integral equation
n(x) = 1 + cj e X I n(x')dx'
0
and evaluating the integral at the endpoints, we obtain two
equations which can be solved for the two unknowns c^ and c^
(3.3) and (3.1)
n(o) = c B
. 1 k
n(o) = 1 + cj e (c cos /~k x + c sin /~k x B ) dx'
0 2 k
(3.3)
(3.1)
From
also
n (1) = c cos/ k + c sin/ k B^
1 2 k
1
J
0
n(1) = 1 + cj e ^ X'( c^cos /k x + c2sin i/~k x B^) dx'
therefore
0 _ Q v I _ A
c B = 1 + c/ e (ccos / k x' + c.sin S k x' B ) dx'
1 1 2
k 0 k
c1cos/r~k+c,.sinVr~it B2 = 1+cj e X ^c^os^k x+c^in/Hc x B2) dx'
k 0 k
Carrying out the integrations in the above equations, we get (cos =
cos ^ k, and sin = sin ^~k for simplicity in writing).
p2 c g ___
c = 1 + c e (Bcos + J k sin) + 3 +
1 k 2Pc
c _ O ^ ___ Q ^ o
e (Psin cos) + ^ k + (e 1)
31
c^os vr~k+c2 sin
__ of c _g
/ k = 1 + c Yfe ( 3cos + /~k sin) Be +
^2 ___  __ Q
Tr ( $sin / k cos) + / k e
2 pc
3 ,, 3,
k ne )
Collecting coefficients of c^ and c2
c 1 ~z e ^ (geos + /~k sin) + 3
1 2 pc
c c 2
e ^ (3sin /~k cos) + i/k = 1 +  + (e ^1 )
c cos  ( Bcos + /k sin) Be +
1 2 p
1 _ / i 3 3 Be B
c sin q (3sin + k cos) + vke = 1 + + (e 1)
2 2p k k
Solving for c^ and c2 we get
C_ =
HHS (X2 + sin + X4)
1 (X2 + sin X4)(1_x4) X2(1X1 cos + X3)
C2 =
C (1X cos + X )
J_______1___________3_
(X2 + sin + X4)
* Q __
where X. = 5 e (Bcos + ^ k sin) + B
1 2 P
1 B
X. = rs e (3sin *r~k cos) + /~k
2 2 P
X.
= i; 3cos + ^~k sin 3e
3
23
X = ^5 3sin ^ k cos + ^ k
4 2 3
_1+Â£t is
k k
3
e
RHS
APPENDIX 2
THE DERIVATION AND SOLUTION TO EQUATIONS (3.5)
From equation (2.4)
e. (x) = e. (x) + aT (n(x)f(x)) +
11 "^ext
(2.4)
(1ct) Â§ fe. ClF1 9 y dx*
2 0 X2 Lxx)2+i3
e. (x) = e. (x) + aT_(n(x)f (x)) +
i_ i . 2
2 ext
(1o) f /e. (x')f
2 o L(x'x)2+i3
dx1
where e. (x) goes to 1 + 2 S2for S=i7 or higher.
W Y+ _ _
Therefore, (2.4) becomes (X = 1a and cos = cos / k, sin = sin / k,
and dropping the subscript i for simplicity)
e (x) = 1 + ff.,1 )S2 + aT1 (n(x)f(x)) + XcJ e w le1 (x1 )dx'
1
e xx'
2 v
__ 8 r 8(xx*)
= R + aT. (c^cos+c.sin 1) + XcJ e e.(x')dx'
112k J 1
^ ) 3(x'x) .., ,
+ J e e (x')dx'
where R = 1 + ^~S2
Using Leibnitz' rule, we differentiate once and get
v 8 (xx')
'e '(x) = ot^ / k (c^sin + c2cos) + n c 3J e e^(x')dx'+
(1 )
33
^(x) + / $e
b(x'x)
e^x' Jdx'e^x)
g(xx')
e ''(x) = otT /~k (c cos + c sin) + Xc3 3_[e
2 1 1 2 0
1 t
(e (x)') + / ge ^e^x')dx'e.,(x)
x
(x')dx' +
(2)
'1
Multiplying equation (2.4) by 3 and subtracting from (2) we get
2 22 q2
3 e (x) e (x) = R3 + aT R (c cos + c sin (1 + )) +
2 2 112 K
(3)
aT^kfc^cos + c2sin) + 2Xc3 e^(x)
We must now insert the expression for e^(x) (from 2.4) into the
above equation and then differentiate two more times. The resulting
expression is
2 iv
3 e *(x) e^ (x) = 2c3ctk (c^cos + c^sinJjT^ + T^X) +
0.2 o ? 3(xx') 3(x'x) '
2(Xcp) 3 J e e^x'ldx' + J e e^x') dx
0 x
2
e2(x)
(4)
After substituting the expression for e^(x) from (2.4) into equation
(3), we must now subtract equation (4) from this, giving us a
fourthordered differential equation for e2(x).
iv 2 4 2 2
e2 (x) 23 e2"(x) + (3 (2 Ac3) Je^x) = (2Xc3)
2c
a (c.cos+c_sin)(T. + TA) aft (T)(T_f3 + 2XcÂ£t )
12 12 2cP 1 2
+ R 32(32+ 2Xc3).
(5)
The righthand side of equation (5) is of the form A cosx + B sinx +
C, where A, B, and C are constants. Therefore, a particular
solution to this, is
e2 (X)
A cos
B sin
k2+34(2*c3)2+ 232k) (k2+34(2*c3)2+232k)
3 (2^c3)
(k is still = 3(2c 3))
2 ____ __
(2c3 )a(T,j+T2^)(c^cos +c2sin)
(32(2c3)2+34 (2Ac3)2 +233(2c3))
a B'
_ (T.3 +2^c3t ) ^ a.2 o2 \ o
2c _J_________2_ R 3 (3 +2Ac3)
2C"3 (34(2Xc3)2) 34(2Xc3)2
After much simplification, this reduces to
(T.+T.^)(c
,, 12 1 2 2c a 1 2 RP
% ~
The complementary solution to equation (5) is
m^x m x m^x m^x
e (x) = be + be + be + be
2 12 3 4
c
where m. = (3  2^c3)
i
1/2
Therefore, e2(x) = e2 + e2 ^x^
c p
4 m.x
= E b.e 1 + Y.cos ^**k x + Y.sin f~k x + Y_
. l 1 2 3
i=1
where Y^ =
(T1 + T2)C1
1 + *
(T1 + T2)C2
y2 rrrrj
23c(1 ^)(3t + 2*cT } r g
3 (32c)(32(2^c)2) P"2Ac
R = 1 + ^i1>s2
The derivation for e^(x) would proceed in the same manner except
that T^ and T2 would be switched.
APPENDIX 3
DERIVATION OF THE CONSTANTS IN EQUATION (3.7)
In order to find the dependence of (eq. 3.6) upon temper
ature, we must substitute the expressions for e^(x) and e^tx) into
the coupled integral equations (eqs. 2.4) and evaluate the integrals
at the endpoints, giving vis four equations and four unknowns. (b^
is related to a^, therefore giving us four unknowns instead of
eight.) As stated in the thesis (p. 17), b^ = a^, b2 = a2, b^ =
a2, b^ = a^, b^ = a^. This relationship can be found by carrying
out the analysis for e^(x) (as was done for e2(x) in Appendix 2) and
comparing a^ with b^. Using the expression for e^x) and e2(x) from
Appendix 2 and evaluating e^(x) at x = 0
1 Bin I
e (0) = A +F (0)+^c/e ' X'e (x')dx* = A +F (0) + ^
I Z ^ I / &
1 a , 1 o , m x'
/e X e^x'Jdx' = A2+F2(0) / e X (a^e + (simi
similar
terms involving a^, a^, and a^)
1 n , 1 fi I 1 O 
+ Je Y^cos / k x'dx' + / e X Y2sin v' k x'd' + / e ^Sf^dx'
where A = R a T (1+)
/, JC
R = i + 2lyDs2
R 1 (Y+1 ) S
F2(x) = a T2(c^cos /"k x+c2sin r\ x)
After doing the integrations and collecting the terms involving a^,
we get (let cos k x = cos; sin x = sin for simplicity in
36
writing)
i[1 * [i   3[i
^ = ^g~ce^(3cos + / ksin) + pj
m 3
Ac(e 1)
m33
m 3
+ j1+Me 1)
4  in 3
Y2 X C / 6 __ I Y A c # _ i
le (3sin / kcos) + vkj ^ le (3cos + ^ ksin) + 3
23c
+ a2 + F2(0)
Now we must evaluate the integral equation (2.4) at x = 1;
e^(1) = A +F (1)+Ac/e ^ ^ X ^(x'Jdx' = A+F(1) +
1
2 2
1
2 2
 _g lx' miX
Ycje I '(a )e dx + (similar terms involving a a ,
0 1 2 3
and a)
4
/e"3l1X,y cos x'dx + / e_f3l1"Xt sin *'r_k x'dx +
0 1 0 2
Again, carrying out the integrations and collecting a^, we get,
[S Ac(e 1e'P)] f
He +=r^J\
m m2 f
2 Ac(e e ,
e t, 1 +
m_
[
m3 + *c(e 3
a31 6 m
rf],[
m4 Xc(e 4aB
3
'f k cos + / k
M1
)>
A2 + F2^1^ + ~2lL(^COS + ^ k sin P e + sin 
3,
(1e )  X^cos X2sin X3<
37
(where , Y^, X^, and X3 are given in text of the
thesis.)
Now, we must evaluate the equation for e^fx) at the same two
endpoints. This will give us two more equations involving, a^,
making a total of four equations and four unknowns.
At this point we could use a matrix solver to solve for a^;
however, as was stated earlier, we would like to find expressions
for the heat transfer which are independent of temperature. There
fore, a^ must be solved by hand so that the dependence upon plate
temperatures, T^ and T^ can be separated out.
Writing out the four equations and collecting terms on the
righthand side that depend on T^ and T^ will give us
(D f + COEF1 ( 1] + a2 Q COEF 1(2 )] + a3 [l COEF 1(3)] +
a4 ^ COEF 1 (4^)j = Y1P1
3
+ Y P + Y_P_ x x +R + a T (c (1 + ))
22 33 13 21 k
(2) a
^ 1 + C0EF2( 1 +a2e 2 C0EF2 (2^j +a3 e 3+COEF2(3)j
0
in
a e + C0EF2(4
3 yi bi
X. + R + T_ ccos + c.sin (1 + )
o 2 1 2 k
+ Y2 W2 + Y3 W3 X1CS X2Sln 
3 2
(3) a1 [l + C0EF1(l7J + a2 Ql COEF 1(2)] + a3 Â£l + COEF loTj +
[l COEF 1 (4 )j = X^t X0P^+ X,P, Y,Y,+ <* T,(c,1 )) + R
2 2 3 3 1 3 1v 1
,m. _ _m
3 2
"k'
ai[e 1 + C0EF2 (1 7j +a2 e 2 C0EF2 (2 )J +a3 Je 3+COEF2(3)j
rm4 n
a e + C0EF2 (4 )J = X1 W1 + x2 w2 + X3 W3 Y cos Y sin 
32
Y + R + a T (c.cos + csin 1 )
3 112 k
(4)
38
t2 + Vr,
where xi = TTT ci
T + XT
X2 = 1 + X C2
X = 2L X c T. + L3T + M + MGS'2
'1
and L =
2 3( 1X)
(32c)(82(2Xc)2)
M =
G =
3
3  2Xc
2(X1)
(X+1 )
and COEF1(1) =
COEF2(1) =
"1.3
Xc(e 1)
111,3
m. 3
Xc(e e )
111,8
COEF1(2) =
m?3
Xc(e 1)
m23
, etc.
m2
C0EF2 (2 ) = ~eQ etc.
m23
Y,, Y2, and Y3 equal X,, X2, and X3 with T, and T2 switched, and P,,
P2' P3' W1' W2' and W3 arS a^'*' *n<3ePen<3ent f temperature and depend
on H only.
After solving equations (1), (2), (3), and (4) for aand collecting
terms that involve T, and one will obtain very complicated
expressions for a^ which can be written symbolically as
(5)
(6)
(7)
(8)
From equations (2.5)
a1 = T1 ri + T2 r2
a2 "5, si + *2 S2 + S3
a3 + ?2 *2
a4 = *1 ui + *2 U2 + U3
q(x) = a (e. (x) T n(x))
1 x2 2
39
m x m x m x m x
= a b^e + l>2e + b3e + b^e + Y cos /""k +
Y sin V k + Y_ T (c, cos (/ k x) + c sin (l/k x)
where ^ = a^ b2 = a2, b3 = a3 b4 = a4# and q 1 = / q^xjdx
0
Doing the integration and writing the final result in terms of the
variables used in equations (5) through (6) we get
(c^sin / k + c2(1cos / k))
_ r mi S m>,
!l = T L Jslli + , + +
a 11 1 m^ 1 m2 1 m3 1 m4
3
[ mi
(e 1
r
2
(1+X) k
1)
32
+ 2Xcl +
m
(e 1) (e J1) (e 41) +
+ s + t + u
(c^sin / k + c2<1cos / k))
(1+X) k
m m
(e 1) (e 1) . 8 . GS2B
f. g ' f q f . r
3 m2 3 m4 32Xc B2Xc
where L was defined previously in this appendix.
The coefficient multiplying T^ is the A(H) (which is a function of H
only) in equation 3.7. The coefficient multiplying T2 is the B(H).
2
The s3 and u3 contain terms with S (S is the speed ratio) and these
are collected together with the
G3
term to give us the C(H) term.
B2Xc
And, finally, the terms left over are grouped as D(H). A computer
program which calculates values of A, B, C, and D is given in
Appendix 4.
APPENDIX 4
COMPUTER PROGRAMS
41
C
C
C
C PROGRAM TO CALCULATE THE FOUR CONSTANTS A,B,C, AND D IN THE EQUATION
C FOR THE FEAT TRANSFER Q.
C
C
PROGRAM GGtOUTFUT, DATA1, DATA2 DAT A3, DAT A4, DATA5 DAT AX,
1TAPE1=DATA1, TAPE2= DATA2, TAPE3= DAT A3, TAPE4= DATA4, TAPE7= DATA5
1 ,TAF8=DATAX)
DIMENSION TEMP(4),TP(4),T22(5),HH(72)
DIMENSION Y(10),X(10),AA(10)
DIFENSION RRM(2,2),AB(2),RA(2),WK(6)
REAL NTRLGYE
PRINT*,' INFUT LAM3A'
READ*, RLAMDA
ALPHA=1R.AFCA
C
C THE SFEED, S, SHOULD BE AT LEAST 7 OVER HIGHER INORDER TO SATISFY
C THE ASSUMPTION OF HIGHT SFEED IN THE DERIVATION OF THE EQUATIONS
C FOR THE HEAT TRANSFER.
C
PRINT*,'INFUT S'
READ*, S
G=.33333333333333
C G = 2*(1GAMMA)/(1 +GAMMA)
PRINT*,' INPUT T1 '
READ*, T1
T22(1)=.1
T22(2)=.5
T22(3)=1
T22(4)=2
T22(5)=5
DO 22 1=0,19
HH(1 + 1) = .01* 1+ .8
WRITE(8,*) HHU+1)
22 CONTINUE
HH(21)=1
WRITE(8,) HHC21)
DO 23 1=1,50
HHC21+I)=1+,1*l
WRITEC8,*) HH(21+1 >
23 CONTINUE
HH(72)=10
WRITE(6,*) HH(72)
DO 999 11=1,5
T2=T22(II)
DO 999 KM=1,72
H=HH(MM)
C PRINT*, H=
C READ*, H
C PRINT*, 'T1=
C READ*. T1
C PRINT*,'T2='
C READ*, T2
C
C BEGIN CALCULATING THE VAULES FOR BETA AND C WHICH ARE USED IN THE
C APPROXIMATION OF THE KERNEL.
C A LEAST SQUARES FCTHOD IS USED.
C
HOLD=0
XSQDSUM=0.
H0LD1=0.
B=0.
C
C CHOOSE NINE POINTS BETWEEN 0 AND 1.
C
42
X(1)=.1
X<2>=.2
X(3) = .3
X(4)=.4
X(5)=.5
X<6)=.6
X(7) = .7
x(8)=.e
X(9) = .9
N=9
DO 20 1=1, N
Y( I )=((H**2)/2)*(1/
20 CONTINUE
DO 2 1=1, N
AA( I )=ALOG(Y (I))
HOLD=AA(l HHOLD
2 CONTINUE
NTRLGVE=HOLD/N
DO 4 1=1,N
H0LD1=XCI J+H0LD1
4 CONTINUE
XB/=H0LD1/N
DO 5 1=1,N
B=AA( I )*X(I )+B
5 CONTINUE
D=NTRLGYE*XB*N
DO 6 1=1,N
XSQDSUM=X(I )**2+XSQDSUM
6 CONTINUE
E=XSQ DSU M (XB AR *2) *N
BZERO=NTRLGVEXB /* (BD)/E
0=EXP(EZERO)
BETA=(BD)/E
C PRINT*,0=',C
C PR INT*,'BETA=i,BETA
BETAfBETA
B=BETA
C
C CALCULATE K WHICH IS BETA TlfS ( 2C BETA)
C
REALK=B*(2*CB)
SQK=SQRT(REALK)
C
C BEGIN CALCULATING C1 AND C2 WHICH ARE USED IN THE SOLUTION OF
C THE NUMBER FLUX EQUAT ION.
C
A1 =EXP(B)*(B*COS(SQK )+SQK*S IN(SQK))
A2=EXP(B)*(B*S I N( SQK)SQK*COS( SQK))
A3=(B*C0$(9QK)+SQK*S IN(SQK))
A4=(B*S I N( 90K)S0K*COS( SQK) >
A5=C/(B**2+REALK)
BB1=1A5*(A1+B)
BB2=A5*(A2+SQK)
BB3=COS(9(?K)A5*(A3EXP(B)*B)
BB4=SIN( SQK)A5* (A4+SQK*EXP (B))
RRM(1,1)=BB1
RRM(1,2)=BB2
RRM(2,1)=BB3
RRM(2,2)=BB4
SUM=1+((B**2)/REALK)+(B*C/REALK)*(EXP(B)1)
AB( 1)=SUM
AB(2)=SUM
M=1
N=2
MA=2
IB=2
IA=2
IJCB=0
oo OOOO ooo o o o
43
LINEQ IS A SYSTEM S LB ROUTINE WHICH SOLVES AN M BY M MATRIX.
CALL LEQ IF(RRM, IA,N,MA,AB, IB,M, IJ0B,WK, IER)
PRINT*,'C1=',ABC 1)
PRINT*, 1 IER=1, IER
PRINT*, 'C2=',AB(2)
C1=AB<1)
C2=AB(2)
NOW THAT BETA,C,C1, AND 02 HAVE BEEN CALCULATED WE CAN BEGIN
CALCULATING THE CONSTANTS FOR THE HEAT TRANSFER.
RM1 = SQ RT (BET A**2+2*RL AM3A *C*B ETA)
RM2= SQRT (BETA**22*RLAMDA*C*8ETA)
RM3=RM1
RM4=RM2
T0P=2*B ETA*C* (1 RL AM)A)
B0T=(BETA2*C)*(BETA**2(2*RLAMDA*C)**2)
RL=TOP/BOT
CHOLD1=( 1EXPCBETA))/BETA
RP3=FLAMDA*C*CH0LD1
D1=C1/(1+RL AM)A)
D2=C2/(1 + RLAMDA)
Q=C1 (1+BETA/ (2*CBETA))
SQK= SORT (2*BETA*OBETA**2)
AHCLD1=(EXP(BETA)*(BETA*COS( SgK)+90K*S INC S?K))
1+BETA)/(2*BETA*C)
RP1=FLAM)A*C*AHOLD1
BHOL D1=(EXP(BETA)* (BETA*S IN(SgK)S0K*COS(SgK))+SgK)/
1 (2*8ETA*C)
RP2=RLAM)A*C*BH0LD1
VAR1=RL*(BETArt2*RLAM)A*C)*(RP31)+(HRLAM)A)*(D1*RP1
1D1+D2*RP2)+ALPHA*0
RM=BETA/ (BETA2*RL AM)A*C)
VAR2=2*G*(RM*(RP31)+1)
VAR3=2*(RM*RP3RMH)
W3=RP3
AHOL D2=(BETA*COS(S0K)+S0K*SIN(SgK)BETA*EXP(BETA))/
1 (2*BETA*C)
W1=RLAM)A*C*AHOLD2
BHaD2=(BETA*S INC SQK)S0K*COS(SQK)+S0K*EXP(BETA))/
1 (2*BETA*C)
W2=RL AM)A*C*BH0LD2
Z=C1*C0S( S0K)+C2*S IN(S0K)( 1+BETA/(2*CBETA))
VAR4=RL*(BETA+2*FLAMDA*C)*(W31)+(1+R.AMDA)*(D1*W1D1*
1 COS(S0K)O2*S INC S0K)+D2*W2)+ALPHA*Z
V AR5=2*G*(RM*(W31)+1)
VAR6=2*(RM*(W31 )+1)
PRINT*, RM1 = ',RM1
PRINT*, RM2=',RM2
TOP 1 = RL AMDA*C* (E XP (RM2B ETA) 1)
BOT1 =RM2BETA
HCLD1=TOP1/BOT1
T0P2=FL AMDA*C* (EXP (RM4BETA)1)
B0T2=RM4BETA
H0LD2=TOP2/B0T2
T0P3=RL AK)A*C*(EXP(RM2)EXP (BETA))
BOT3=RM2+BETA
H0LD3=TDP3/B0T3
TO P4=RL AM)A *C* (EXP (RM4)EXP (BETA>)
44
B0T4=RM4+BETA
H0LD4=TOP4/B0T4
TOP5=R.AM3A*C*(EXP(RM1 BETAJ1)
B0T5=RM1BETA
HOL D5=TO P5/ BOT5
T0P8=RLAH)A*C*(EXP(RM3)EXP(BETA))
B0TB=RM3+BETA
H0LD6=T0P8/B0TB
T0P6=RL AM)AC*(EXP (RM3BETAJ1)
B0T6=RM3BETA
HOL D6=7D P6/ BOT6
T0P7=RL AM)AC* (EXP (RK1)EXP (BETA))
B0T7=RM1+BETA
HOLD7=TDP7/BOT7
RL5=2*( 1+HCLD5)
R_6=2*(1+H0LD6)
RL7=2* (EXP(RM1 )+H0LD7)
RL8=2*(EXP(RM3)+H0LD8)
DENB=RL5*RL8RL6*RL7
COF1 =R_ (BETA2*RL AM)A*C )* (RP3+1)+(1 RL AW)A )* (D1 *RP1+
1D1+D2*RP2)ALPHA*0
C0F2=RL*(BETA2*RLAFDA*C)*(W3+1)+(1RLAFCA)*(D1*W1+D2*
1W2+D1 *COS( SQK )+D2*S IN(SQK) )ALPHA*Z
RL1=2*(1H0LD1)
RL2=2*( 1H0LD2)
FL3=2* (EXP (RM2) HOLD3)
RL4=2*(EXP(RM4)H0LD4)
HOL DA= (VAR 1 *RL 4 RL2* V AR4)
HOL DB=(S**2)*V AR2+VAR3
HOL DO (S**2) *VAR5+VAR6
DENA=RL1*RL4RL2*RL3
A2=(T1*HOL DA+T2*H0L DA+RL4*H0L DBRL2*H0L DC) /DENA
R1 = (OOF 1 *RL8 RL6*00F2)/DEI'B
R2=R1
S1 =( VAR1*RL4RL2*YAR4) /DENA
S2=S1
SMALLT1 =(C0F2*RL5RL7*C0F1) /DENB
SMALLT2=SMALLT1
U1=(RL1*VAR4RL3*VAR1)/DENA
U2=U1
U3=(R_1*(VAR6+VAR5*S**2)RL3*(VAR3+YAR2*S**2))/DENA
TP( 1 )=RM1
TP(2)=RM2
TP(3)=RM3
TP(4)=RM4
DO 150 1=1,4
TEMPI I ) = (EXP(TP( I) )1 )/TP( I)
150 CONTINUE
RNUM=C1*SIN(9}K)+C2*(1C0S(SQK))
DN=(1+RLAMDA)*SQK
RNDN=RNUM/DN
CON ST A= R1 *TEMP (1)+S1 *TEMP (2)+SMAL LT1 *TEMP (3)+U1 *TEMP (4)
1RNDN(BETA/(BETA2*C))+R.AMDA*2*C*RL
CONSTB=R2*TEMP(1)+S2*TEMP(2)+SMALLTZ*TEMP(3)+U2*TEMP
1(4)+RNDNRL*BETA
FRSTPAR=(RL4*VAR2RL2*VAR5)*TEMP(2)/DENA
SCNDPAR=TEMP(4)*(RL1*VAR5RL3*VAR2)/DENA
TH RDPAR=G*B ETA/ (BETA2*RL AM3A*C)
CONST OFRST PARf S ON DPARf THRDPAR
FRTHPAR= (RL4*VAR3RL2*VAR6)*TEMP( 2) /DENA
FTH PAR= (RL 1 *V AR6 RL3 V AR3) *TEMP (4) / D EN A
SXTHPAR=BETA/(BETA2*RLANEA*C)
CONSTD=FRTH PAR+ FTH PAR+SXTH PAR
C PRINT*,'00NSTA=',CONSTA
C PR I NT*,1 CONSTB= *,CONSTB
C PRINT*, OONSTO',CONSTC
C PR I NT*,1 CONSTD= CONSTD
o o o o
45
Q1 =ALPHA*(T1 *C0NSTA*T2*C0NSTB+( S**2)*CONSTC<CONSTD)
PRINT*,'T1=',T1
PRINT*, T2=',T2
PRINT*,'lt=', H
PRINT*, 'Q1 = ',Q1
IF(I I.EQ.1) WRITEd,*) 01
IF(II.EQ.2) WRITE(2,*) Q1
IF(I I.EQ.3) WRITE(3,*) Q1
IF(II.EQ.4) WRITE(4,*) Q1
IF( I I.EQ.3) WRITE(7,*) 01
999 CONTINUE
END
oooo ooooo oooo ooo
PROGRAM TO CALCULATE THE L I FT ON PLATE 1.
PROGRAM GGG(0UTPUT,DATA1,DATA2,DATA3,DATA4,
1DATA5, TAPE1=DATA1, TAPE2=DATA2, TAPE3= DATA3, TAPE4= DATA4,
1TAPE7=DATA5)
DIMENSION X(7),WT(7)
DIMENSION HH(20) ,TT(5)
PRINT*,' INPUT H'
READ*,H
PRINT*,' INFUT T1,T2'
READ*, T1 ,T2
Pf=3.1415927
TT(1)=.1
TT(2)=.5
TT(3)=1
TT(4)=2
TT(5)=5
HH(1)=.8
HH( 2)=1
DO 111 1=1,16
HH(l+2>=1 + l/2.0
11 CONTINUE
THE INTEGRAL EQUATION THAT CALCULATES THE LIFT IS EVALUATED USING
GAUSSIAN QUADRATURE.
WCX) VALUES #?E THE WEIGHT FUNCTIONS.
X(1 >=.0254460
X(2)=.129234407
X (3 > = .29707 7424
X(4)=.5
X(5) = .702922575
X(6)=.6707 65592
X(7) =.974553 956
WTC1>=.0647424831
WT(2> = .13985 26 957
WT(3)=.1909150253
WT(4)=.208979591 8
WT(5)=. 1909150253
WT(6)=. 13985 26 957
WT(7)=. 0647424831
00EF=(H**2)*(2.0/3.0)*(1/PI)
PRINT*,' INPUT T1 '
READ*,T1
THE SFEED, S, SHOULD BE AT LEAST 7 OR HIGHER. SEE PROGRAM FOR
CALCULATION OF HEAT TRANSFER FOR EXPLANATION.
PRINT*,' INPUT S'
READ*, S
PRINT*, INPUT LAM3A'
READ*, RLAMDA
DO 999 II 1=1,5
T2=TT(III)
DO 999 MMW=1,20
H=HH
C0EF=(H**2)*(2.0/3.)*(1/PI)
SUM=0.
DO 1 K=1,7
U=X (K)/SQRT( (X (K )**2) +H**2)
UPRIME=( 1X(K))/SQRT((1X(K) >**2+H**2>
PARTI =SQRT(TR (X (K) ,2, H, T1,72, RL AKOA, S))
PART2= SQRT (TR(X (K>, 1, H, T1, T2, RLAMDA, S))
VAR1 =PART1 PART2
TRM1 =2*GAMMA(1X(K),H)/( H**2)
TRM2=1 / GAMMA (1 X (K), H)
TRM3=EL PTK (UPR I >E) / ((1 X (K))*G AMMA (1 X (K), H))
TRM4=2*GAMMA (X (K), H) / (H**2)
47
TRM5=1/GAMMA(X(K) ,H)
TRM6= EL FIX (U) / (X (K) *GAMMA (X (K), H))
TRM7=PI/(H**2)
Y=R'IPiaX(X(K) ,H)*VAR1*( (1/ (1XCK)) )*(TRM1 TRM2)*
iaPTE(UHRIPE)TRM3+(1/X(K))*(TRM4TRM5)*aPTE(U)TRM6TRM7)
C WHERE RNHFLX(X) IS THE VALUE OF NUP6ER FLUX AT POINT X.
C a PTE IS THE EL I FT I CAL INTEGRAL OF THE FIRST KIND.
C ELFTK IS THE a I FT I CAL INTEGRAL OF THE SECOND KIND.
C TR(X) IS THE TEMPERATURE AT POINT X, WHICH WILL BE A
C FUNCTION OF THE ENERGY AT X AT THE NUPBER aUX AT X.
C
SUM=WT(K)*Y+SUM
1 CONTINUE
C
C CALCULATE THE LI FT ,
C
a IFT=OOEF*SUM
C PRINT*,1 a IFT=*,a I FT
C PRINT*, 'T2=',T2
C PRINT*, >H=',H
IF(I I I.EQ.1) WRITEU,*) FLIFT
!F( I I I.EQ.2) WRITE(2,*) a I FT
IF(I I I.EQ.3) WRITE(3,*) FLIFT
IF C11 I.EQ.4) WRITE<4,*) a I FT
I F( I I I.EQ.5) WRITEC7,*) FL I FT
999 CONTINUE
END
FUNCTION aFTK(X)
C
C ELPTK(X) IS THE aiPTICAL INTEGRAL OF FIRST KIND.
C
DIPENSION Y(0s50)
N=10
Pl=3.14159
B=PI/2
AfO
DX=(BA)/(2*N)
DO 1 K=0,2*N
Y(K)=1/(1(X**2)*(SIN(A+K*DX)**2))
Y (K)=SQRT (Y (K))
1 CONTINUE
S=(DX/3)(Y(0)+4*Y(1))
DO 2 1=2,N
S=SHDX/3)*C2*Y(2*l2)+4*Y (2*11))
2 CONTIMJE
S=S+(DX/3)*Y (2*N)
apix=s
RETURN
END
FUNCT ION aFTE(X)
c
c ELIFTICAL INTEGRAL OF THE SECOND KIND.
C
DIPENSION Y(0:50)
N=10
Pl=3.14159
CO
rM
s.
a. o
co
cn
ft?.
W CN
*L:
>oÂ£tuj
^ CN
3 _L
? ft
? ^
O &
* *ro
CN ~
z
a
?
CO
S
X CO
VI* LU ^ CN N U4 >
QJto U X 3 X CO
r~2Â£9ol*Â£*EP
>&S&S&dÂ£
CN
CN
I
+
X
*
X
X
d
X
A
X
X
d
ii
CO
<
d
p
CO
<
E
a!
_2
i*
a*
x
<
E
d
X
CN
X z x
w LU
1
si!
e <.
Z I
2 V
bS
ozÂ£
o
^a:
*3385
f r
O
O 5 G3 u5 ^*iu
IZCNZ
LU t LU Lii
q: ^
eg 3C lu *
m
Â£*&Â£Â£Â£
I
m 2. LU ^ UJ ^
X w  yj  CO Q I Q
Â§Â£5Â£a!.Â£Â£uJ{EGjSÂ£3
O
o
X
d
X
Â§
z
o
X
X
CN <
CN
X
Zi
d
s
. a' o
~ ^ CN
NT *
^ S.
X CN
n *
w ^ CN
Â£ S3
>:>
o ;
00 c
:s
I to ;
2 fÂ£ ^ ^
z 2 Â£Â£Â£a'2"l
q q 5 aÂ£ tb
co
co
co
co
co
CO
CO
CO
K\
CO
CO
CO
CO
O CO
I?
o a
x
d
x
z
X
ooo
< o
1 Â£
< O =3
TASq ..........
^ O (N co O psffl 0>
^xxxmxxxxxxxxx
II II II II II II u u u
X
u
CD
o
o o o o o
20
N=9
DO 20 1=1,N
Y(I)=((H**2)/2)*(1/(X(I)**2+H**2))**(1.5)
CONTINUE
DO 2 1=1, N
AA(I)=AL0G(Y(I)1
HOLD=AA(l 1+HOLD
2 CONTINUE
NTRLGVE=H0LD/N
DO 4 1=1,N
H0LD1=X( I1+H0LD1
4 CONTINUE
XBAR=H0LD1/N
DO 5 1=1,N
B=AA( I )*X( I )+B
5 CONTINUE
D=NTRLGVE*XBAR*N
DO 6 1=1,N
XSQDSUW=X( I )**2+XSQDSUM
6 CONTINUE
E=XSQDSUMCXBAR*2)*N
BZER0=NTRLGYEXB/*(B0)/E
OEXPCBZERO)
BETA=(BD)/E
PRINT*,'0',C
PRINT*,'BETA=',BETA
BETA=BETA
B=BETA
REALK=B*(2*CB)
SJK=SQRT(REALK)
A1 =EXP(B1* CB*C0S( SJK)+SQK*S INC SpK) 1
A2=EXP(B)*(B*SINC 9pK)SQK*COS( 9JK))
A3=(BCOS( SQK)+SQK*S INC SpK))
A4=(B*SIN(SJK)S?K*C0S(SQK))
A5=C/(B**2+REALK)
BB1 =1A5* (A1+B)
BB2=A5*(A2+9QK)
BB3=C0S( SQK1A5* (A3EXPCB) *B)
BB4= SIN(S?K)A5* (A4+SJK*EXP (B))
RRMC1,1)=BB1
RRM(1,2)=BB2
RRMC2,1)=BB3
RRM(2,2)=BB4
SUI4= 1 +(C B**2) /REALK)+( B*C/REALK>* (EXP (B) 1)
ABC 1)=SUM
AB(2)=SUM
M=1
N=2
HA=2
IB=2
IA=2
IJ0B=0
CALL LEQIFCRRM, I A, N, MA, AB, IB,M, IJCB,WK, IER)
FAINT*,'C1 = ',AB( 11
PRINT*, IER= 1, IER
PRINT*,C2=,AB(2)
C1=AB(11
C2=AB(2)
AAA=BETA/(BETA2*C)
RNMFLX=C1*C0S( SQK*
RETURN
BID
oo o OOO
FUNCTION ENERGY(QQ, I IPAR, H, FF, CC, ALPHA, S)
FUNCTION TO CALCULATE THE ENERGY AT FOINT QQ.
DIMENSION TEMP(4),TP(4)
DIMENSION Y(10),X(10),AA(10)
Dl KENS ION RRM(2,2),AB(2),RA(2),WK(6)
REAL NTRLGVE
IF(I IFAR.EQ.I) THEN
T1 =FF
T2=CC
ELSE I F( 11 FAR. B?.2) THEN
T2=FF
T1 =CC
ELSE
ENDIF
RLAKOA=1ALPHA
G= .33333333333333
G = 2*(1GAMMA)/(1+GAMMA)
HGLD=0
XSQDSUM=0.
HCLD1=0.
B=0.
X(1)=.1
X(2)=.2
X(3)=.3
X<4)=.4
X(5)=.5
X(6) = .6
X(7)=.7
X(8) = .8
X(9)= .9
N=9
DO 20 1=1 ,N
Y(I)=((H**2)/2)*(1/(X(I)**2+H**2))**(1.5)
>0 CONTINUE
DO 2 1=1, N
AA( I )=ALOG(Y(I))
HOLD=AA(l J+HOLD
! CONT IMJE
NTRLGVE=HOLD/N
DO 4 1=1,N
HOLD1=X(IJ+HOLD1
I CONTINUE
XBAR=HOLD1/N
DO 5 1=1,N
B=AA( I )*X( I
> CONTINUE
D=NTRLGVE*XBAR*N
DO 6 1=1,N
XS(JDSUM=X(I )*2+XSQDSUM
5 CONT IMJE
E=XSQDSUM(XBAR**2)*N
BZERCNNTRLGV EXB AR* (BD)/E
O=EXP(BZER0)
BETA=(BD)/E
PRINT*,'O',C
PRINT*,'BETA=',BETA
BETAfBETA
&=BETA
REALK=B*(2*CB)
SQK=SQRT(REALK)
A1=EXP(B)(B*C0S(SQK)+SpK*SIN(SQK))
A2=EXP(B)*(B*S I N( SQK)S0*COS( St?K))
A3=(B*COS( SQK)+SQK*S IN(SQK))
A4=(B*S I N( SJK)SQK*CO$< 9QK))
A5=C/(B**2+REALK)
BB1=1A5*(A1+B)
51
BB2=A5*(A2+9QK)
BB3=C0S( SQK)A5*(A3EXP(B)*B)
BB4= SIN < SQK)A5* (A4+SQK*E XP (B))
RRM(1 1)=BB1
RRM(1,2)=BB2
RRM(2,1)=BB3
RRM(2,2)=BB4
SU M= 1 +((B**2) /REALK)+(B*C/REAL K)* (EXP(B) 1)
AB( 1 )*=SUM
A8(2)=SUM
C CALL GAUSS(3,2, RRM, AB)
M=1
N=2
MA=2
IB=2
IA=2
IJ0B=0
CALL LEQ IF(RRM, I A, N, MA, AB, IB,M, IJCB,WK, IER)
C PRINT*,'C1=',AB(1)
C PRINT*, IER=', IER
C PRINT*,'C2=',AB(2)
C1=AB( 1)
C2=AB(2)
RM1 = SQ RT (BETA *2+2 *RL A>CA *C *B ETA)
RM2= SQRT (BETA**22*RLAMDA*C*8ETA)
RM3=>RM1
RM4=RM2
S=10
T0P=2*B ETA*C* (1 RL AMDA)
B0T=(BETA2*C)*(BETA**2(2*RLAMDA*C)**2)
RL=TOP/BOT
CH0LD1=( lEXP(BETA) )/BETA
RP3=RLAMD1 *C*CH0LD1
D1=C1/(1+FL AM)A)
D2=C2/(1+BLAMDA)
Q=C1( 1+BETA/ (2*CBETA))
9QK= SQRT (2*BETA*CBETA**2)
AHCLD1=(EXP(BETA)*(BETA*COS(SQK)+SQK*SIN(SQK)J+BETA)/
1 (2*8ETA*C)
RP1=RLAM)A*C*AH0LD1
BH0LD1=(EXP(BETA)*(BETA*S IN( 9QK)SQK*C0S( SQK) )+SQK)/
1 (2*8ETA*C)
RP2=BH0LD1*RLAM)A*C
V Afi 1 =RL (BETM2*RL AfdA*C) (RP31)+(1+RL AM)A) (D1 *RP1
1D1+D2*RP2)+ALPHA*Q
RM=BETA/ (BETA2*RL AK)A*C)
X1=D1*T2+R.AMDA*01*T1
X2=D2*T2+RL AM3A *0 2*T2
X3=2*RL *RL AM)A*C*T 1+RL *B ETA*T2+RMt RM*G (S**2 j
V AR2=2*G* (RM* (RP31)+1)
VAR3=2*(RM*RP3RM+1)
W3=P3
AHOLD2=(BETA*COS( 9QK)+SQK*S IN( SQK)BETA*EXP(BETA))/
1(2*BETA*C)
W1=FLAFDAC*AH0LD2
BH0LD2=(BETA*S IN(SQK)SQK*COS( SQK)+SQK*EXP(BETA>)/
1(2*BETA*C)
W2=R_ AM)A*C*BH0LD2
Z=C1*C0S( 5QK)+C2*S IN( SQKM 1+BETA/(2*CBETA))
V AR4=RL*(BETA+2*RLAMDA*C)*( W31)+(1+PLAMDA)* (D1*W101*
1 COS (SQK)02*S IN (SQK >+D2*W2)+AL PH A*Z
V AR5=2*G*(RM*(W31)+1)
VAR6=2*(RM*(W31)+1)
C PRINT*, 'RM1 = ,RM1
C PRINT*, RK2=',RM2
TOP1=RLAMDA*C* (EXP(RM2BETA)1)
52
B0T1 =RM2BETA
H0LD1=TOP1/B0T1
T0P2=R.AM)A*C*(EXP(RM4BETA)1)
B0T2=RM4BETA
H0LD2=TOP2/B0T2
T0P3=RL AM)A *C*(EXP (RM2)EXP (BETA))
B0T3=RM2+BETA
HOL D3=TO P3/ B0T3
T0P4=RLAW)A*C* (EXP (RM4) EXP (BETA))
BOT4=RM4+BETA
H0LD4=TOP4/B0T4
T0P5=RLAM)A*C*(EXP(RM1BETA)1)
B0T5=RM1BETA
HOL D5=TO P5/ B0T5
TOP8=RLAM)A*C*(EXP(RM3)EXP(BETA))
B0T8=RM3+BETA
H0LD6=TOP8/B0TB
T0P6=RL AFDA *C* (EXP (RM3BETA) 1)
BOT6=RM3BETA
HOL D6=TD P6/ B0T6
TOP7=RL AM)A*C*(EXP (RM1) EXP (BETA))
B0T7=RM1+BETA
HOL D7 =TO P7/ BOT7
RL5=2*( 1+H0LD5)
RL6=2* (1+H0LD6)
RL7=2* (EXP(RM1 )+HCLD7)
RL8=2* (EXP (RM3)+HOL D6)
DENB=RL5*RL8RL6*RL7
COF1 =FL (BETA2*RL AKDA*C)(RP3+1) +(1 RL AK)A ) (D1 *RP1 +
1D1+D2*RP2)ALPHA*Q
C0F2=RL (BETA2*RL AM)A*C)(W3+1)+(1 RL AH3A) (D1 *W 1+D2*
1W2+D1 *COS( SpK)+D2*S I N( SQK))ALPHA*Z
RL1 =2* (1HOLD1)
RL2=2*(1HCLD2)
FL3=2* (EXP (RKZ) HOL D3)
RL4=2*(EXP(RM4)H0LD4)
Ha DA= (VAR 1 *RL4 RL2* V AR4)
Ha DB=( S**2) *V AR2+V AR3
Ha DC= (S**2) *V AR5+V/K6
DENA=RL1*RL4RL2*RL3
A2=(T1*HaDA+T2*HaDA+RL4*HaCBRL2*HaDC)/DENA
R1 = (COF 1 *RL8RL6*C0F2) /DEhB
R2=R1
S1 =( VAR1 *RL4RL2*VAR4) /DENA
S2=S1
SMALLT1 =( COF2*RL5RL7*COF1)/DENB
SMALLT2=SMALLT1
U1=(RL1*VAR4RL3*VAR1)/DENA
U2= U1
TP(1)=RM1
TP(2)=RM2
TP(3)=RM3
TP(4)=RM4
DO 150 1=1,4
TEMPO )=(EXP(TP(I ))1)/TP(l)
150 CONTINUE
A1 =R1*T1+R2*T2
S3=( HaDBRL4RL2*HaDC)/DENA
A2= S1 *T1+S2*T2+ S3
A3=SMALLT1T1+SMALLT2*T2
U3=(HaDC*RL1HaCB*RL3)/(DENA)
A4=T1*U1+T2*U2+U3
ENERG Y= A1 *EXP (RM1 0g)+A2*EXP (RM2*QQ)+A3*EXP (RM3*QQ)+A4
1 EXP(RM4*QQ) +X1 C0S( SQK*QQ) +X2*S I N( SQK*QQ) +X3
C
C
RETURN
END
