Bamberger, et al. 5. This approach essentially ...... Heinrich,. J. C., and. Zienkiewicz,. O. C.,. "Natural convection in a Square. Enclosure by a. Finite-element,.
IP
•/F _t" /
.
NASA
Technical
Memorandum
r
103721
A Finite Element Model of Conduction, Convection, and Phase Change Near a Solid/Melt Interface
Larry Lewis
A. Viterna Research Center
Cleveland,
January
Ohio
1991
j.-
#
A
FINITE
ELEMENT PHASE
MODEL
CHANGE
OF
NEAR Larry
National
CONDUCTION, A
CONVECTION,
SOLID/MELT
A.
AND
INTERFACE
Viterna
Aeronautics and Space Administration Lewis Research Center Cleveland, Ohio 44135
ABSTRACT
Detailed flow
is
These a
understanding
required
for
systems
often
of
accelerations
range
of
many
heat
transfer
aerospace
include
thermal
phase or
and
systems.
change
effective
fluid
and
operate
over
gravitational
fields. An which
approach
requires
conservation as the
an
the laws
equation
A
is
program.
The to
variables,
presented Galerkin
solve
the
along
with
algorithm.
used
the
velocity.
a
procedure
marching for
are
Newtonian
fluid.
spatial
implicit
energy
Lagrangian
solving
in finite of
the
elements
two are
well
form
of
governing a
computer
element the
method
field
Crank-Nicolson Lagrangian
and
as
two-dimensional
the
variation
Quadratic
internal Linear
in
the
the mass,
developed
of
presented
property
implemented
form
an
and
variable
for
and
is of
momentum, The
for
systems
solution
energy,
state.
coordinates numerical
such
simultaneous
equations
equations
used
analyzing
of
of
governing
Cartesian
is
to
time
elements
are
components
of
used
for
the
pressure. The location the
of the
temperatures
are determined
internal
energy
general
in that
phase change,
it
This
can describe
heat
phase change with
Analytical
calculated
approach
is quite
transfer
without
a sharp
interface,
from this
and
researchers
model are compared to
studying
transient
convection,
and phase change and are found
agreement.
The numerical
significant
computer
procedure
resources,
when compared to similar Several
as well
an interface.
results
of other
interface
from the
and pressure.
phase change without
those
solid/liquid
methods are suggested
by other
to reduce
times.
ii
to be in good
presented
but this
studies
conduction,
is
requires not
unusual
researchers. the
computational
as
ACKNOWLEDGMENTS
I would Professor help to
in
Gene
E.
completing
acknowledge
doctoral
especially
the
Smith and
like and
to Dr.
critiquing
assistance
of
committee.
iii
thank Robert this all
the
co-chairmen
Siegel
for
work.
I also
the
members
their
of
wish my
TABLE
ACKNOWLEDGEMENTS NOMENCLATURE
OF
CONTENTS
...................
iii
......................
v
CHAPTER I.
INTRODUCTION
2.
LITERATURE
3.
PROBLEM
................. REVIEW
4
............
12
.............
25
Statement
Governing 4.
...............
FORMULATION
Problem
1
NUMERICAL
Equations
APPROACH
Subdividing the Domain Formulating the Element Equations Assembling the System Equations Solving for the Transient Response 5.
RESULTS Case Case Case Case
AND i: 2: 3: 4:
Computer 6.
CONCLUSIONS
VERIFICATION
..........
72
Conduction only Phase Change by Conduction Buoyancy-Driven Convection Combined Conduction, Convection, and Phase Change Resource Usage AND
RECOMMENDATIONS
.......
105
APPENDICES
......................
109
REFERENCES
......................
172
iv
NOMENCLATURE
a C D e Fo g
i J Ja k K L n m M N Nu n P Pr q r Ra T t u v V W x Y
P o T
thermal diffusivity specific heat (also capacitance geometric dimension internal energy Fourier number acceleration (also node
gravity) number
time Jakob
step increment number
matrix)
thermal conductivity stiffness matrix latent heat unit normal vector unit mass
tangential matrix
vector
interpolation Nusselt number number
of
(shape)
global
functions
nodes
pressure Prantl number heat flux number of nodes Rayleigh number temperature time
per
element
velocity in x-direction velocity in y-direction velocity Gauss-Legendre weighting coordinate in cartesian coordinate in cartesian coefficient of thermal transformed coordinate
factors system system expansion
absolute viscosity density normal stress shear stress (also dimensionless time) kinematic viscosity field variable field value at a node dimensionless transformed
location coordinate
of
v
phase
interface
e
v
parameter in time-marching recursion (also dimensionless temperature) prescribed nodal value divergence of a vector
Subscripts
1 S 0
Matrix
[J { } [ ]
and
classifies classifies classifies classifies classifies classifies
superscripts an area or volume integral a surface integral an approximate value of field liquid state solid state initial condition or reference
notation single single matrix
algorithm
row matrix column matrix
vi
variable
state
CHAPTER
1
INTRODUCTION
The
need
heat
transfer
more
complex
true
in
systems
fluid
aerospace
must
operate
or
industry
and
effective
conditions
of
problems,
change
the
time
most
ground-based
of
are required
low-gravity
analytical
over
a
For
systems
range
of
systems
under
expensive.
phase
change
are
not
and
must
these
numerical
no
fields.
many
studies
or
Such
and
facilities
or
little
fluid/thermal
for
laboratories.
particularly
available.
difficult
designs
reliable
where
operata
modeling man
is
highly
gravitational
experimental
Earth-orbiting
is
as
This
where
data
for
increasing
environments
investigations
gravity
is
tools
devices.
in
phase
Experimental
predictive
flow
experimental
include
Because
and
analytical
thermal/fluid
accelerations
low
general
the
supporting often
for
possible be
done
reasons
method
in on
a
would
be
very
valuable. The purpose analyzing undergoing
intent
of
numerical the
heat
phase
this
research
approach transfer change.
and and Such
is
to
computer fluid an
develop
a general
program flow
analytical
of
for materials
tool
would
2
significantly and
aid
in
Many tool
reduce our
such
systems
and
casting
methods.
analysis in
lunar
the
space
As
would
working
latent for
and
are
well
developed
recently,
computer
known
been and
on
might
in
Klann
the
this
heat would
dark
periods
low
in
device
the to
be
Aeronautics Freedom
or
space
power
energy
to
engine.
A
solar
cycle
device
is
Station
34,
have
provide of
the
heat
energy Space
Burnsl0). and
fluid
and
to
numerical
problems for
modeling
heat empirical
geometries.
been
years. transfer
easy
for
techniques
have
many heat
generally
programs
convective
simple
they
thermal
National
Space
Brayton
during
for
limited
on
concentrate
development
programs
mostly
a
and
computer
effect
by
approaches
thermal
research
its
the
storage
(see
Analytical
of
of
system
orbit,
modeling
system
and
thermal
power
Station
storage
and
and
application
(NASA)
collect
as
a thermal
described
fluid
heat
the
power
management
processing
freezing
Another
for
fluid
batteries
examined
results.
a computational
material
inadvertent
Administration
base.
system
be
space.
such
cryogenic
as
required
experimental
of
such
could
experiments
the
advanced
designs
the
Space
modeling
or
of
of
of
Systems
initial
temperatures
and
as
devices
withstand
of
applications
analyzing
management
number
understanding
practical
exist,
used
the
to
modeling
transfer,
the
subject
Today, by
apply. fluid however,
relationships Today,
for
computational
conduction Until flow
and
had based
on
to
3
convection
using
numerical
differences
and
handle
complex
more
commercially and
in
additional
restricts
to
Chapter with
phase
become
available
most
are
of
these
research
applications.
widely
varying
properties
present
methods
work
is
change
needed
to
fluid/thermal
are
The
phase
most
to
oriented
having
of
analyze
thesis
the
presented
phenomenon
for
develop
such
general
problems
and in
recommendations
for
discussed
and
to
Chapter
presented
of
5.
The
in
which
further
work
PHASTRAN,
3 and
equations
verification
6,
of
developed in
the
the
thermal,
Chapters
governing
Chapter
program,
presented
solutions
with
computer
ends
is
problems.
Results
are
review
numerical
of
approach. model
2 a
change
development
is
exhibit
finite
with
change.
dealing
and
of
as
Some
particular
further
methods
have
however to
such
geometries.
application
Hence
In
and
scope
that
the
problems.
phase
flow
complexities
materials
purpose
elements
available,
limited
and
finite
methods
literature
fluid
flow,
4 cover
the
and
the
the
approach
main
numerical
body
concluding are
appendices.
of
the
remarks
given.
during
and
this
A research
CHAPTER
2
LITERATURE
Much
research
materials
undergoing
association
with
examples.
and More
modeling
58
change
on
related
and
Luther,
including
the
Thornton
method; and
To related to
numerical
discussed The
32,
classical
in
and
this
Wilkes
are
the
flow.
the II
modeling
of by
convective
61
on
thermodynamics.
of
the
analytical
those
methods
those
change
to
4,
Huebner
element transfer;
important
to
the
the
to:
problem;
heat
most
related
4
referred
finite
first
approach
is
Baker
on
phase
detailed
numerical
the
for
a
change
2
the
in
literature
reader
method;
topic,
important
interest
For
phase
on
are
the
on
some
food,
environment.
in
the
of
its
an
65
Larsen
Sonntag
of
been
space
difference
research
followed
the
on
Zienkiewicz and
has
fluid
38,
analysis
industries
there
and
and
summarize to
The
subjects
finite
Arpaci
VanWylen
applications.
research
Lunardini
Carnahan,
and
of
the
because
processes
phase
discussion
to
semi-conductor
examples
modeling
devoted change
recently,
these
Many
been
phase many
metallurgical,
Stefan
has
REVIEW
works
papers
related
problem
will
fluid
flow.
phase
change
be
5
problem
will
not
mathematical complex
discussed.
derivation, geometry
properties
is
classical
and
and
of
involved
finite
the
finite
superiority
rather
the most
finite
element Otis
region
of
into
solved
the
finite
space
assumed
uniform
and
latent
works
melting
the
include
to is
have
a
lesser
not
because
method, of
have
the
been
within heat
intervals. each
of
but
two
methods.
devoted
was
required
in
of
a
of
materials
initially
and
41
pseudo
time
by
to
element
modeled a
dividing
Temperature
volume
effect
method
problem
The
analyses
on
change
and
difference
source. terms
phase
This
development recent
papers
the
approach.
46
the
the
of
material
change
method,
finite
its
54.
of
method.
chronological
Indeed,
Savino
difference
the
phase
models
element
of
of
in problems
noteworthy
and
numerical
to
dependent
Some
Siegel
elegant
application
modeling 9,
the
the
extent,
its
impractical.
Kreith
Most
Though
temperature
analytical
Budhia
any
be
at
as
a
coordinate
variable
and or
the was
any
instant
moving
heat
transformation was
limited
finally
at
to
the
fusion
temperature. Murray the
interface
location
differential The
front
at
the
location
suggested
was
equation
differential
balance
Landis
for
equation phase was
calculated
approach by
the
velocity
of
was
derived
from
front. then
an
set
The equal
solving the an
temperature to
the
by
which a
phase
front.
energy at
freezing
the
new
6
temperature. Springer approach Again the
to
the
remainder
the
the
and
involved
interface.
and
fluid
of
need
Only
researchers
convection
transfer. field
Such
procedures
solved included
difference
formulation and
equations. track
compared
stream
Again the
in
phase
favorably
the of
the
Murray
front. with
track
the
best
method
phase in
Chapter the
experimental
3.
effects
of
heat in
the
solution
problem used
conservation
Tien's
for
times.
He
form
the
change.
radiative
iterative
and
finite
nonlinearities
fluid.
a
the
balance,
included
change
function
the
the
computational phase
to
formulation
follows
or
from
a
energy
was
fluid
equal
the
as in
explicitly
have
require
the
the
introduce
increased
convection
vorticity
to
in
enthalpy
their
this
the
which
and
Tien 59
in
set
conduction.
temperature
method
effects
equations
heat
conduction
on
few
was
determined
of
multidimensional
a
was
Landis
dimensions.
front
used
maintained
this
two
53
to
and
temperatures
for
approach
the
They
in
Because
discussion
natural
to
the
Sparrow
integral
Murray
phase
and
instead
eliminated
analysis
the
formulation. an
the
front
solution
variable
difference
method
at
solid
difference
used
phase
temperature of
dependent
56
temperature
Shamsundar
More
Olson
track
fusion
finite
and
of
a
using
a
was
used
momentum
approach
numerical data
natural
finite
laws
the
Landis
with
results on
the
freezing
7
of
naphthalene. Valle
60
also
included
solution
but
solved
method.
The
conservation
the
stream
effects in
function and
terms
of
an
of flow, his
that
interface
motion
of
the
melt
the
interface
more
and
the
limitation
Brown
of
disappearance,
not
implicitly
the date
seem
the
most
detailed
and
included
effects.
and
to
the
track based
and
18
have
shape
an
and
being
able or
finite
field
approach
variables However,
mesh to
problem
boundaries,
elegent
techniques. moving
the
deforming
the
fixed
is
other
merging
and
on
transformation.
transformed
This
approached
easily
fragmentary
O'Neill
I
used
a
method
so
that
of
which
which allowing
this
formulations,
has
handle distribution
phases. Albert
Valle
however,
approaches
have
moving
regions one.
with
as
coordinate
solution
as
not
works
using
interface
efficient
approach,
several
solid is
did
heat
at
Tien,
of
Landis.
and/or and
of
effectively
problem
Ettouney
to
terms
latent
radiation
work
approach
as
of
problem
the
in
fluxes
one
and
element
formulated
heat
his
finite
The
was
in
developed
were the
change
to
recently,
grids
couples
of
tension,
and
change
element
phase
this
Murray
More phase
motion
This
results
concluded
that
front
surface
were
the
temperature.
imbalance
the
convection using
laws
interface.
analyses
compared
problem
and
phase
solid/liquid
fluid
the
natural
of
transfinite
of
the
8 mappings
in
finite
conjunction
element
technique,
phase
front
was
Again
there
is
Ettouney
and
have
studied 52
finite enthalpy the
of
algorithm If
adjusts
the
one
iteration
is
significantly
also
affect
of
fluid
the used
25
the
finite difficult boundary
used
to
vary
of
problems years.
finite
flow layer
method
one
computational
numerical
flows
including and
even
flows. has
been
materials
to
variable More
developed
modeling
progress
the
years,
the
to
viscous
many
flows,
property
handle
over
solutions
recently, to
with
the
solutions
slow
but
interface.
Over
provided
only
times
computer-based
has
problems
spacing,
remarkable
method.
on
nonlinearities.
methods
made
the
iterations.
for the
of
Depending
grid
the
difference
(thermo-hydrodynamic) element
within
near
the
Schneider's
the
has
method
53"
convergence
especially
using
variation
Sparrow
of
Initially,
difference
a
converge
rapidly
application flow
last
moves
accuracy,
that
The
only
researchers
problem
interface,
number
reduces
the
properties
the
application.
schemes.
with
and
the
associated
some
change
along
of
the
interface
This
phase
Shamsundar
movement
costs
intensive
the
for
its
problem,
the
approach.
above
computational
technique of
mesh
restricts
change
formulated
method
amount
high
mesh
tracking
fixed
mentioned
numerically
difference
boundary-moving in
a
which
phase
less
to
limitation
the
moving
improvement
method
the
a
compared
the
of
modeling
Schneider
made
Brown
Because with
with
flows the
many
finite of
the
may
9
same to
problems. be
The
particularly
geometries
and
Gallagher,
et
application
boundary al. 21
finite
prevented
the
problems.
use
Later,
methods
broadened
variety
of
fluid
Olson
44
Galerkin
3
applied to
originally
currently
the element
Hood formulated
and the
finite
used
of
the
problems. element
The
method
Navier-Stokes
of
equations, to
weighted
the
application
of
finite
to
exact however,
practical
of
to
methods
lack
application
flow residual
elements
to
a
pseudo-variational
the
formulation
weighted
to
developed
a in
used
(discretized) 31
technique
flow.
nondiscretized
widely
Taylor
residual
incompressible a
most
The
approach
method
of
of
Galerkin is
formulating
the
equations.
also
used
Navier-Stokes
the
Galerkin
equations
formulation;
the
approach
function.
Comparison that
and
variational
elements
shown
complex
3,
the
formulation;
formulation. suggests
of
a
velocity/pressure vorticity
complex
equations.
the
viscous
criteria,
finite
Baker
finite
incompressible stream
22
the
been
involving
examples
to
often
of
applied
the
Baker
of
has
problems.
two-dimensional of
many
element
forms
problems
element
problems
variational
terms
contain
applications
the
in
method
conditions.
finite
continuum
derive
element
useful
of
Early some
finite
and of
the the
three
stream
purely
these
velocity/pressure
in
and
ways:
function
stream
three
criteria
the
and
function
formulations
formulation
may
have
10
several
advantages.
dimensions. stress
Pressure,
boundary
appears
It
to
is
readily
velocity,
conditions
require
less
extended velocity
can
be
three
gradient,
easily
computational
to
and
handled.
time
And
than
the
it
other
formulations. Recently,
more
considerations these
for
nonlinear
solution
been
good
problems
Important
of
has
obtaining
fluids
conditions. choice
attention
quality
over
aspects
of
technique,
given
to
solutions
a wider
this
element
the to
range
include
of
proper
types,
and
mesh
refinement. Gartling, equations
et
in
terms
the
convergence
two
element
flow
properties
those
system
applicable
convergent.
between
an
8-node
direction
reduced
adequate of
the
solution
most
mesh rapid
was
and
studied algorithms,
used
Laminar
to
represent
techniques,
they
the
full
unsymmetric
and
more
generally
counterparts. procedure
significant
with
element
refinements.
walls
quadrilateral
element, its
mesh
symmetric
No
pressure
plane
superior
finite
severalsolution
Newton-Raphson
rapidly
of
of
solved
their
the
quadrilateral
Of
the and
several
which were
than
particular,
Finally,
and
problem.
that
because
velocity
converging
a nonlinear
equation
of
types,
between
found
al.23"formulated
the
In
was
the
difference
element
and
8-node
being
was a
preferred
in
formulation
refinement
was
required
of
the
found
13-node
complexity
variation
most
solution
and in
use.
the field.
ii
Ben-Sabar choice
of
ratio
of
found
that
and
boundary
boundary of
Reynolds
implicit
use
of
conditions
20
the
the
effect
problems
terms
are
the
velocity
are
necessary
instabilities
developed
an
element
terms
viscous
of
where large.
and
with
the
the They
surface
to
dominate
with
equivalent
the
finite
for
and
flow
an
alternating
method
compressible
comparison
delay
the
increasing
element
flows
applied
past
a
direction where
the
method
rectangular
finite
approach
the to
object.
difference
to
be
In
scheme,
he
computationally
efficient. Solutions be
the
problems
in the
Though
beyond directed
systems
of
methods
will
computer processing
to
subject
require
are
on
diffusive
numerical
finite
convection
to
investigated
number.
Fletcher
more
to
consistent
appearance
6
conditions
convective
traction
found
Caswell
coupled
fluid/thermal
of
research,
much
three-dimensional expenditure the at
scope improved
nonlinear be
architectures capabilities.
continue
particularly
space.
Such
transient
problems
often
of
significant
computer
resources.
of
this
a
of
solution
equations.
dependent
problems
on
study, methods The
the
providing
number to
use
of
availability vector
and
efforts
solving many of
of new
parallel
large these
CHAPTER 3
PROBLEM
FORMULATION
Problem
The
problem
transient
selected
temperatures,
velocities,
and undergoing
in
and
solid
satisfying contained shown
in
in
force and
or
3.1.
fluid
nor
is
crystallization motion
is
heat
of
due
reference
a
frame
It such
exist
fluid
is as
could
and
to
can
model.
conditions
induced
eutectic
properties
geometry
flux,
or
material
state
arbitrary
heat
include velocities.
gravitational and
is
both
body inertial
included.
energy
(surface
transfer
and
mechanisms
are
also
restricted
The
the
fluid
substance
variable
Boundary
are
free
supercooling
fluid
is
effects
surface
included, of
of
accelerating
viscous No
with
determine
rates,
pure
change.
temperatures, the
a
equation
vessel
Figure
prescribed Flow
a
analytically
transfer
in
states
general
in
to
heat
phase
fluid
a
is
pressures
material
statement
not
tension)
by
radiation.
of
nucleation
considered.
to
12
laminar
The
are
effect
or In
flow.
effects
addition,
the
13
Insulated
_:iiiii
Pr:aStribedFlux
iiiii_ i!iiiiii_ii! iiiii!_iii_i_:_i_i_i
:::._i!
Figure
3.1
as
LaErangian
reference
frame
equations, In particles space.
one
Example
the
science
or
Eulerian adopted. these
Lagrangian
which The
in
of
can
Phase
and
independent
Change
engineering
depending To two
System
Equations
on
formulate approaches
approach, be
TP;eS_rbtdur e
}iiii::
Governing
Problems
_
all
can
be
classified
the
viewpoint
the
governing
must matter
identified
as
they
variables
in
the
be
adopted.
consist move
or
of
through
Lagrangian
system
14
are
x0,
which
Y0' a
In by
z0
and
specified the
continua
t
where
fluid
approach,
of
quantities.
field
are
the
To
derive
the
governing
on
one
the
in
governing
control
laws
volume
equations. in
and
approach
thermal
adopted
The
solving
physical
laws
as
well
is
problems, undergoing To
and
our
a
time.
attention
volume.
with
If
we
apply
differential differential
which
analysis
are
formulated
modeling
the
phase
equations
most
problems
and
change
expressing
of
energy
2.
Conservation
of
momentum
3.
Conservation
of
mass
a
thermodynamic this
expressed Such
characterized
z
governing
Conservation
equation.
heat.
to the
particularly
commonly
approach
y,
focus
to
i.
as
energy
x, we
of
to .
independent
control
set
time
is
the
the
problem three
of:
Because it
problem
at
are
The
the
coordinates
through
processes
a
a
the
here.
solution
includes
passed
called
the
are
coordinates
obtain is
z0
equations,
of
we This
fluid
spatial
space
Y0,
element
Eulerian
variables
area
x0,
problem
phase further
is
important The in
formulations however,
equation
they change discuss
of
dominated to
of
are may at
a
this,
of
thermal the
energy
temperature
quite be
by
consider
conservation terms
state.
valid
discrete let
us
form is
and for
inappropriate
aspects, the
most specific
single for
temperatures. consider
of
two
phase
materials
15
situations, the
one
other
will
in
in
which
which
a sharp
interface
a mushy
region
with
a
sharp
interface
top
and
is
no
formed
sharp
and
interface
exist. An
example
of
water
with
in
Figure
of its
3.2.
Suppose
bottom
the
might
sides
water
was
be
a thin
insulated
layer
as
initially
at
shown a
Initial Temperature,
To=m f
CIL)/_EC , ÷(~CHKEF)/DEC . ÷(~CHKTF)/_EC , ÷0 . _EC:TIME÷TCTL[1] , INITIAL , NTS÷I . CIC÷IOpO a TCTL[3]÷TCTL[3]e2 n _ND:MSG+'TIME STEP CHANGED STATUS MSG TCTL[3] . ÷0 . _TOP:STATUS 'NOT CONVERGED' ÷ .
ON
WITH
.
OR
O'S
MAKE SURE TABLE IS A MATRIX CONVERT ROW TO A MATRIX FIND LARGEST OF COLUMNS RESHAPE TABLE ADD ROW(S)
ABILITX
TO'
BLANKS
.
TO
CONVERGE
AND
FIELD
VALUES
EXIT IF TIME STEP CAN'T CHANGE CHECK FOR TOO MANY ITERATIONS IF ENERGY OUT OF RANGE. DECREASE IF TEMP. OUT OF RANGE. DECREASE EXIT. TIME STEP IS OK SET BACK TIME TO BEGINNING REINITIALIZE PROBLEM DECREASE TIME STEP COUNTER REINITIALIZE ITEM. COUNTER DIVIDE TIME STEP B_ 2 MESSAGE TO USER DISPLAY AND RECORD MESSAGE EXIT MESSAGE TO USER END EXECUTION
127
AGAIN [0] [1] [2] [3] [5] [63 [7] [e] [g] [zo] [11]
AINTCON [0] [I] [2] [3] [_] [5] [6] [7] [8] [g] [10] [11] [12] [13] [15] [16] [17] [ZS] [ig] [20] [2Z] [22] [23] [2_]
NIT+AGAIN ,
ITER
INTERATIVE
SOLUTION
NIT+ITER+ITER+I
(FOR
NON-LINEARITIES)
,
+(ITER.CIL)/O SETSTAR ,
ITERATION CHECK IF SET LAST
,
ENERGy FLOW UPPROP +(CHKCNV SETSOLID NIT+AGAIN
, NIT)/O A ITER
,
CHECK CONVERGENCE SET VELOCITIES
n
RECURSIVE
NI+I+ANI+N ETA_(NZ.NI)pNI$(Q![I;;])[ANI|] AXEI÷(2.(NIxNI))O(._ETA).(.ETA) LAS÷LSHAPE AXEI . QAS+QLGSHP AXEI A QJAC÷RZE JACOB QAS QAIF+MDET QJAC JACCHK QAIF 81÷RZE[|_I 3 5 7]BF B2+RZE BF QAS , 2
1
,
2
D÷_ 1 2 DXT÷DXT.cD D+_ I 2 DYT+DYT.cD
TIME
STEP
ITERATIONS ITERATION
VARIABLES
ENERGY EQUATION TRE FLUID EQUATIONS BASED ON PF AND OF TF. UF. SOLID NODES
OF
EF
VF AND PF TO ZERO
CALL
AINTGRT ORDER OF GAUSS-LEGENDRE ETA COORDS. XI AND ETA COORDS. LINEAR SHAPE FUNCTION QUAD. LAGRANGIAN SHAPE FNS. JAC08IAN FOR QUAD. ELEMENTS Q_AD. AREA INTEGRATE FACTOR
,
LAS
,
3_(1.NE.oNSHP)pNSHP÷LAS[1;;],
NSHP+W 2 I 3_(I.NE.pNSHP)pNSHP÷QAS[1;;], NST÷NST.cNSHP . D+_ I 2 3_(I,pD)pD+BI[;;I;] . DXT+cD D÷4 I DZT+cD
EACH
COUNTERS TO0 MANY CONVERGENCE
FORM AND SOLVE FORM AND SOLVE UPDATE PROPERTIES
AINTCON N;BI;B2_D;NI;ETA;QJAC;NSHP;LAS:QAS INTEGRATION CONSTANTS FOR USE WITH
NSHP÷_ NST÷cNSBP
AT
3_(1.pD)pD÷BI[;;2:] A 3_(1.pD)pD+B2[;;1;] , 3_(1.pD)pD+82[;;2;] ,
.
CHECK LINEAR QUAD. SHAPE PUT IN SHAPE ADD TO DERIV. PUT IN DERIV. ADD TO DERIV. ADD 20 DERIV. ADD
TO
JACOBIAN GRADIANT MATRIX GRADIANT MATRIX FNS. ELMNT. TYPE 1 GLOBAL VARIABLE FNS. ELMNT. TYPE 2 GLOBAL VARIABLE WRT. X ELMNT. TYPE GLOBAL VARIABLE WRT. Y ELMNT. TYPE GLOBAL VARIABLE WRT. X ELMNT. TYPE GLOBAL VARIABLE WRT. Y ELMNT. TYPE GLOBAL
VARIABLE
I 1 2 2
128
AINTCRT [0] I÷AINTGRT A;WE;WEX;N1 [1] . INTEGRATES FUNCTION [2] . [3] NI÷ANI+I . [.] [5] [63 [7] [8] [9] [10] [113 [12] [13]
AREA tO] El] [23 [3] [_3
BDY t0] [I] [2] [3] [5] [6] [73 [8] IS] [103 [11] [12] [133 [1_] [15] [16] [17] [18] [19] [2o] [21] [22] [23] [2_] [25] C26] [27]
OVER
AREA
÷(lep,A)/_CALAR . ÷(^/(pQAIF):2+OA)/_XT R STATUS EEE_GEz] , ÷0 . _XT:A÷Ax(20tppA)W(2_pA)pQAIF ÷QALC . _CALAR:A+AxQAIF QALC:WE÷(N1.N1)pNI+(_[2;:])[ANI;] WEX÷(oA)p_((x/(I+pA)).(I*pA))p(.WE)x°_WE I÷+/[I]WEXxA ,
AREA;RN n CALCULATES
AREAS
SAREA+GFxSINTGRT VOL÷AINTGRT 1 n
AND
OF
IN
, ,
VOLUMES
OF
FROM
COOR.
SYSTEM
ELEMENTS SIDE AREAS VOLUME OF
R÷HDY FT;B;M;NC;N;D;S . EVALUATES R VECTOR
XI-ETA
1+ORDER OF INTEGRATION CHECK IF A IS SCALAR CHECK IF pA IS LIKE pQAIF MESSAGE TO USER EXIT MULT. 87 INTORT. FACTOR JUMP TO FINISH INTEGRATION MULT. 8Z INTGRT. FACTOR WEIGHTING FACTORS RESHAPE WEIGHT FACTORS NUMERICAL INTEGRATION
.
1 ,
8÷FT FSTCM 8C , NC*8[;I] R+(_.(SNI+I).NE)pO
ELEMENT
BOUNDARY
AND
OF ELEMENTS ELEMENTS
INTERNAL
CONDITIONS
BC'S FOR THIS FIELD LIST OF BC TYPES INITIALIZE R VECTOR
,
TYPE TO
ZERO
R
&00P:÷(O=pNC)/NXT1 N÷I+NC m NC+(NzNC)/NC +(N=l)/&00P . B÷N FSTCM 8 . D+RN FSTCM 8 . ÷(O:xlpD)/_XT1 . ÷(N:2)/_2
A
.
M÷'BOUNDARY CONDITION STATUS M . ÷&OOP . _2:R÷D 8DY2 R . +&OOP .
NOT
DEFINED
(BDY)'.
LOOP ON BOUNDARY CONDITIONS TAKE FIRST 8C TYPE REMOVE THIS TYPE FROM LIST IGNORE PRESCRIBED 8C'S BC'S FOR THIS BC TYPE BC'S FOR THIS REGION NUMBER EXIT IF NO SUCH BC'S 80UNDARY CONDITION TYPE 2 WARNING TO USER DISPLAY AND RECORD MESSAGE CONTINUE LOOPING PRESCRIHED FLUX 8C END OF LOOP
R
NXTI:÷(1 2=FELT[FT])/&INEAR._UAD STOP R _UAD:S÷GSS[1;;] . +NXT2 . _INEAR:S÷LSS[1;|] . NXT2:R+. 2 3 I_((-I÷pS),pR)pR R÷3 1 2 _W((pR)DS)xR .
n
n
CHECK IF LINEAR ELEMENTS WRONG ELEMENT TYPE QUADRATIC SHAPE FUNCTIONS FINISH CALC. OF R VECTOR LINEAR SHAPE FUNCTIONS RESHAPE R AND MULTIPL_ BY SHAPE FNS.
129
8DY2 [0] [1] [2] [3] [.] [5] [6] [7] [8] [g]
G÷D 8D_2 R;DI , MODIFIES ELEMENT R MATRIX FOR BOUNDARY COND. TYPE 2 (PRESCRIBED FLUX) A D[;1] IS THE SIDE NUMBERS , D[|2] IS THE PRESCRIBED FLUX VALUES A ÷(O:I÷DD)/0 . EXIT IF NO 80UNDAR_ CONDITIONS G+R . INITIALIZE RETURN VARIABLE DI+c[2]D . NESTED ARRAY OF BOUNDARY COND. 8DY26SUB"D1 . HANDLE EACH BOUNDARY CONDITION G÷R n RETURN R VECTOR FROM 8D_SU8
8DY2_SU8 [0] 8DYT_SU8 DI=EL:S|PV [1] n SUB-FUNCTION OF BDY2 TO PRESCRIBE EACH (ELEMENT. SIDE COMBINATION) [2] . D1 IS A NESTED ARRAY OF PRESCRIBED SIDE BOUNDARY CONDITIONS [3] A R IS THE R MATRIX FROM 8DY2 WHICH IS MODIFIED BY THIS FUNCTION [5] [6] [7] [8]
BF [0] [1] [23 [3] [_]
CHKCNV [0] [I] [2] [33 [5] [6] [7] [8] [9] [10] [113 [12] [13]
S÷DI[I] n EL÷RSELMT[S:] . PV÷((pR)[2].pEL)_"o"'DI[2] R[S;;EL]÷PV _
_
8÷RZ 8F SICIRS , CALCULATES THE FIELD VARIABLE GRADIANT RS÷I*pS[2:_] 8÷2 I 3 _(NE.pC)pC÷S[2;|].[I.5]S[3;;] B÷(INV RZ JACOB S)INPROD B
SIDE AFFECTED ELEMENTS AFFECTED PRSCRBD. VALUES AT MODIFY THE R MATRIX
INTERPOLATION
INTG.
PTS.
EL
MATRICIES
CNV÷CHKCNV NIT " CHECKS CONVERGENCE OF EF. U_. VF AND PF WITB LAST ITERATION VALUES . CNV RETURNS 1 IF ALL CONVERGED 0 IF NOT ALL CONVERGED 4 CNV÷ERR CBKCNV_SI 'EF' , CHECK EF JUMP IF NO FLOW CALCS. +(FC=O)I_ND . CNV÷CNV^ERR CRKCNV_SI 'UF' . CHECK UF CNV÷CNV^ERR CRKCNV6SI 'VF' . CHECK VF CNV÷CNV^ERR CHKCNVaSI 'PF' A CHECK PF END:+(CNV:O)/&IMIT , CHECK FOR CONVERGENCE STATUS 'CONVERGED WITB' NIT 'ITERATIONS' , MESSAGE TO USER ÷0 A EXIT &IMIT:_(NITI_pTMP)/O RGB+MINMAXt"$"TMP[;I] A RGB+(MEAN RC8)+-I lx(l+TOL)xO.5x-/@RG8 RGT÷MINMAX TF RTC+^/RGTopX)/YL n MIL+~(_ppX)eAS +(~^/(pY):(_MIL)/pX)/ERR1 RT+(MIL/_DpX),(_MIL)/,ppX Y_RT_((pX)[RT])pY ÷APP YL:MIL÷~(_ppY)eAS n ÷(~^/(pX)=(~MIL)/pY)/ERR1 RT÷(MIL/IpoY).(~MIL)/IppY X÷RT_((pY)[RT])pX n APP:R+X FN Y n ÷0 A ERRI:OES 'RANK ERROR'
PHASTRAN [0] PHASTRAN A [1] n MAIN CONTROL FUNCTION [2] . [33 FRONTEND A n [_] TIMESTEP _ _ [5] 8E_ULT_÷RESULTS n [6] STATUS TSTOP n
ELEMENT T_PES 1-ENERGY)
OPERATOR
.
n n
FOR
SELECT
PHASE
ALONG
FROM
GLOBAL
NESTED
ARRAY
AXES
RIGHT ARGUMENT CONTAINS AXES AND Y IF RANKS ARE SAME APPLY THE OPERATOR FIND LARGEST RANK MISSING INDEX LOCATIONS CHECK FOR RANK ERROR RESHAPE AND TRANSPOSE INFO. RESHAPE Y JUMP TO APPLY THE OPERATOR MISSING INDEX LOCATIONS CHECK FOR RANK ERROR RESHAPE AND TRANSPOSE INFO. RESHAPE X PERFORM THE OPERATION EXIT
CHANGE
ANALYSIS
MODEL
UPFRONT. ONCE ONLY FUNCTIONS TIME INCREMENT FUNCTION FORMAT FINAL RESULTS DISPLAY COMPUTER USAGE
PROPDATA6SEL [0] R÷PROPDATA6SEL N [1] n RETURNS SUBSET OF PRQPDATA BASED ON STATE CONDITIONS AND [2] _ N IS NUMERIC VECTOR CORRESPONDING TO SF NOTATION [3] . [_] R+_BQ_Z_ . ALL PROPERTY DATA [5] R+(v/R[3;;]eN)/[2]R n ONLY REQUESTED STATES
PROPERTIES
155
PRSCRB [o] [1]
PRSCRB;A;BNP;KB
[2] [3] [_] [5] [6]
n n
PRESCRIBES PFEQ[1;] PPEQ[2;] PFEQ[3;]
FINITE ELEMENT E@UATIONS POSITIONS IN KBAR (1THRU CORRESPONDING PRESCRIBED (0 FOR SOLID NODES, 1 FOR
IN THE I+pKBAR) VALUES BOUNDARY
PFEQ÷PFEQ[;APFEQ[I;]] n PFEQ+(RDUPL PFEQ[1;])/PFE@
[7] [8] [9] [lO] [11] £129
A÷_I÷pKBAR . PFEQ+(PFEQ[1;]eA)/PFEQ BNP÷~AePFEQ[I;] A A÷-(-BNP)/KBAR n XB+((pA)_PFEQ[2;])xA KB++/PFEQ[3;]x[2]KB RBAR÷(BNP/RBAR)+BNP/KB KBAR+BNP/[1]BNP/KBAR
[13] [I_] [15]
PRSIDE [0]' [I] [2] [3] [4]
B÷BC[;3 2 B÷RN FSTCM
/5] [6] [7] [8] [g]
I
u B
5]
[10]
n
DEFINES
PRESCRIBED
A
AND
ORDER (LEAVE 1ST) IN KBAR
BOUNDARY
CONDTIONS
o A
-DUPLICATE NODE BC'S -LEAVING ONLY THE FIRST RESET GLOBAL VARIABLE PFEQ
n
PRSIDEAADJ [0] [1] [2] [3] [_]
n n
AN÷ON PRSIDEAADJ ADJUSTS ORIGINAL
TYPE;PETS NODE POSITIONS
PFTS+(-I++/^kTYPEzCFT)÷FT AN+ON++/NNG PFTS.
PRSIDEAFLD [0] B PRSIDE6FLD [1] . PRESCRIBES [2] FSTCM .
TYPE_N;V BOUNDARY
[3] [_]
8+TYPE +(O:oB)/O
8
[5] [5] [7] [8]
N+TYPE SIDENODES 8[;1] V÷,W(@pN)pi"i"8[:2] (N V)+((,N)V)PRSIDE_TEMP N+(,N)PRSIDEAADJ TYPE
[9]
PFEQ÷PFEQ,N,[I]V,[0.5]I
IN
A
A
MATRIX
PREVIOUS _DJUSTED
FOR
CONDITIONS
n
A
FOR FIELD NUMBERS
FIELD
TYPE
TYPES (POSITIONS)
FIELD BOUNDARY CONDITIONS FOR EXIT IF NO BC'S FOR THIS FIND NODE NUMBERS FIND ASSOCIATED VALUES
A TYPE n
RBAR
CONDITIONS)
REORDER BC (REG. TYPE FT SIDE ONLY THIS REGION ONLY PRESCRIBED BC'S PRESCRIBE BC'S FOR EACH FIELD MODIFY PFEG TO ELIMINATE
n
PFEQ+DI,[I]D2
KBAR
REMOVE ANY OUT OF RANGE B00L. POS. NOT PRSCRBD. COLUMNS OF K'S PRESCRIBED g'S x PRESCRIBED V_LUES ZERO SOLID NODES NEW R VECTOR NEW K MATRIX
A
B÷I FSTCM B a ((pFT)ocB)PRSIDEAFLD"FT END:DI+PFEQ[I:] DI+¢(C÷RDUPL DI)/DI+_DI D2+¢C/¢PFEQ[2 3;]
VARS.
PUT IN ASCENDING REMOVE DUPLS. ALL POSITIONS n
PRSIDE FT;B;C;D1;D2 o MODIFIES PPEQ WHICH A
GLOBAL
.
HANDLE TEMP. ADJUST MATRIX ADD PREVIOUS
TYPE TYPE
BC'S DIFFERENTLY POSITIONING PRESCRIBED VALUES
VAL.)
TYPE
156
PRSIDE6TEMP [0] D+NF PRSIDEATEMP TYPE|I_C_I1;I2:N_V CONVERTS TEMPERATURE BOUNDARY CONDITIONS TO ENERGY BC'S CZ] [2] , [3] D÷NV . INITIALIZE RETURN VARIABLE [_] ÷(TYPE_I)/O , EXIT IF NOT TEMPERATURE BC NODES AND VALUES [5] (N V)÷NV , SETUP RETURN VARIABLE [6] D÷N.[O.S](pN)pO , [7] II+IEA1 I÷(1 FVG8 SF)[N]el _ 6 . IND. OF NODES IN SINGLE PHASE [83 I÷(PROPDATA6SEL 1 _ 6)[_ 1 2_;] . SINGLE-PHASE PROPERT_ DATA USE TEMPERATURE AND PRESSURE [g] C+V[II].[I.5](1 FVGB PF)[N[II]] , INTERPOLATE SINGLE-PHASE POINTS [10] D[2;I13÷(C INTQR I)[1|1:] n RETURN NESTED ARRAY [11] D÷c[2]D
PRSPRES [o] PRSPRES [1] . PRESCRIBES PRESSURE [23 . [3] ÷(2:QNC 'PRNODE')/NP [_3 [5] PRSP1 [6] ÷0 [7] NP:PRSP2
NODEM A
PRSP1 [0] PRSPI:A_8_NN [1] . PRESCRIBES PRESSURE NODEM: [2] . USES PRESSURE VALUE FROM [3] . [4] NN÷NNG 2 . [5] A÷NNp((I+2xND)pl 0).(ND+I)p0 [6] B÷v/(2.NN)o(NN+z2xNN)_PFEQ[I;] [7] A_(AzO)/A+(~AxB)xA\,(ND+I)*2 [8] 8÷PF[I÷A] A [g] PFEQ+PFEQ.3 lp((I÷A)+SwNN).B.I
PRSP2 [03 PRSP2_N [1] . PRESCRIBES PRESSURE NODEM: [23 . [33 N÷(+/NNG 2 3)+PRNODE[1] . [_3 PFEQ÷PFEQ.3 lpN.PRNODE[2].I
CHECK
IF
PRESCRIBE
PRESSURE PRES.
PRESCRIBED
CASE LAST
_ n .
(TOTAL
PRESSURE
WHERE TOTAL ITERATION
.
NODE
PRESCRIBED MASS
NODE
MASS
CONSTRAINED)
(OPEN
SYSTEM)
CONSTRAINED
NUMBER OF QUAD. NODES (GLOBAL BASIS) BOOL. WHERE PRES. NODES OCCUR BOOL. WHERE VELOCITIES PRESCRIBED PRES. NODES WHERE VEL. NOT PRSCRBED. USE 1ST UNPRSCRBED. PRES. NODE PRESCRIBE THE PRESSURE NODE
CASE
OF
OPEN
SYSTEM
,
MATRIX LOCATION OF NODE NUMBER PRESCRIBE THE PRESSURE NODE
157
PRSVLC [03 CE+PRSVLC FT [1] A RETURNS LOCATIONS OF [23 A [33 CE÷(+/NNG FT)pl . [_3 CE[(PFEQ[1;])]÷O .
QLGSHP [0] [I] [23 [33 [4] [5] [6] [7] [8] [g] [Z0] [11] [12] [13]
PRESCRIBED
SRP+QLGSHP XE;A;C:ETA;XI CALCS. QUADRATIC LAGRANGIAN
VALUES
FOR
VECTOR INSERT
SHAPE
FUNCTIONS
XI÷XE[I:] ETA÷XE[2;] n C÷(ETAxXI.2),(XIxETA.2),[I.5](XIxETA).2 C÷I.XI,ETA,(XIxETA),(XI*2),(ETA.2),C A÷Q_&_xOPAX 2(0 A÷A[:2 5 _ 7 I 8 SHP+SHP,[O.5]C+.x_A A+C_&_xOPAX 2(0 A÷A[:3 . 6 8 7 1 SHP÷SHP,[I]C+.x_A
I 0 I 2 0 2 1 2) 3 g 63 . 0 I I 0 2 I 2 2) g 2 5] n
A
.
QUAD_PLI [0] P+QUAD_PLY XI;X;Y [1] A FORMS POLYNOMIAL FOR 2 DIMENSIONAL QUADRATIC [2] . XI IS A 2 COLUMN MATRIX OF X I VALUES [3] . [4] X÷XI[:1] [5] Y÷XY[:2] n [6] P÷I,X,(X*2),Z,(Y.2),(XxY),(yxX.2),[1.5](XxY.2)
THESE
FIELD
TYPES
OF 1'S FOR ALL NODES O'S WHERE PRESCRIBED
AND
DERIVATIVES
XI COORDINATES ETA COORDINATES BUILD MATRIX OF .COEFFICIENTS CALCULATE SHAPE COEF. FOR DERIV. REORDER DERIV. WRT XI COEF. FOR DERIV. REORDER DERIV. WRT ETA
POLYNOMIAL FUNCTIONS WRT XI
WRT
X VALUES Y VALUES FORM POLYNOMIAL
QUADRGXY [o] Z÷D QUADRGX_ XI;C:K;DX;DY;PXY [1] PERFORMS QUADRATIC REGRESSION ANALYSIS IN 3 DIMENSIONS X Y Z [2] . XY IS A 2 COLUMN MATRIX OF X AND Y VALUES TO BE INTERPOLATED [3] . D IS A 3 COLUMN MATRIX OF X Y Z DATA USED IN THE INTERPOLATION [_] . Z RETURNS INTERP. VALUES AND THEIR DERIVATIVES p (l+pXI) 3 [5] A [6] K+QUAD6PLY D[|I 2] POLYNOMIAL OF INTERP. DATA [7] K+DE:3]|K REGRESSION ON DATA [8] PX)[+QUAD_PLI XI A POLY. OF XY'S TO BE INTERP. [g] Z+PXY+. xK n INTERPOLATED VALUES [zo] C+1 2 0 1 1 2 0 0xK[2 3 1 6 8 7 1 13 COEF. FOR DERIV. WRT. X [11] DX÷PXY+. xC DERIV. OF Z WRT. X AT XY [12] Z+Z, [1.5]DX . CATENATE TO RETURN VARIA8LE [13] C+1 1 1 2 0 2 0 0xK[_ 6 7 5 1 8 1 1] COEF. FOR DERIV. WRT. Y [1_] DY÷PXY+. xC . DERIV. OF Z WRT. Y AT XY [15] Z+Z,DI . CATENATE TO RETURN VARIA8LE
ETA
158
RCOND [0] [1]
_
RC+RCOND RETURNS
SUB-VECTOR
FOR
[2] [3] [_]
RC÷+/[3]((DX RC+RCxOPAX(1
[5] [6] [7]
RC++/RCxOPAX(I 2 _)((DX RC_((ORC),I)oRC A RC_-NODEV AINTGRT RC
RDUPL [0] [11 [2] [3]
R÷RDUPL . n
n
A
BOOLEAN
FOR
X÷IND REDUCE REDUCES (BY
1)
TF)
n
.
DERIVATIVES MULT. BY MULT. BY RESHAPE INTGRT.
REDUCING
DUPLICATE
OF TEMP. CONDUCTIVITY DERIV. SHP.
FNS.
TO COLUMN VECTOR GLOBAL NODE BASIS
VALUES
LEAVING
ONLY
THE
1ST
WITH
MULTIPLE
INDICIES
REORDER INDICIES SEQUENTIALLY AND THE CORRESPONDING VALUES LOCATIONS OF DUPLICATES JUMP IF DUPLICATE INDICIES EXIST RETURN VARIABLE EXIT
A n
C+VAL+I#VALx(-I+pVAL)_C3*I C+((-I+DVAL)¢C3)/C X÷D REDUCE C n
REMOVES
VECTOR
n
[12] [133 [1_1
REMDUPEL [0] R+REMDUPEL
A
n
C3÷C2+(~C2)x((pC2)-I)_C2 D÷C3/INDxC3 n
[3]
FVEB
n
X÷VAL,[O.5]IND +0 n NXT:
n
I),DY
VAL;C;D;C2;CS;I SUMMATION)
[I0] [111
[1] [2]
3)(1
l_
