 

3"; h «ﬁshaizﬁnaiv

 

 

..W.,,H..

.HV...
vr.-

{‘11
-‘uno' 4.4:.

£-..o r.
.. m;

'_“,n». .izi.
.m,‘
.3

‘" "3;“! M ”‘3. ~
‘1

.‘.-,..;-."-;.« 34,.
.....,....;,.

 

mt‘mumnm...

. ..., - u-.. hbﬁﬁu .1.

 

 

' "va‘

ICHIGAN STA

It itit\\\\\\itvtii\itiml

3 1293 00794 3412

 

 

 

 

 

This is to certify that the

dissertation entitled
Crossflow Permeation of Viscous and
Viscoelastic Liquids Through Arrays of

Circular Cylinders
presented by

Craig A. Chmielewski

has been accepted towards fulﬁllment
of the requirements for

 

 

 

 

Ph . D . degree in CHE
11/18/91

Date

 

MS U is an Afﬁrmative Action Hiqual Opportunity Institution 0- 12771

 

 

 

LIBRARY
Mlchlgan State
UntversIty

r... ‘- _-'

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE

 

 

I e

 

 

 

 

 

 

 

 

 

 

 

 

 

‘ﬂﬁﬁ‘l
F__Jl__JL_
_._r—I:

MSU Is An Affirmative Action/Equal Opportunity Institution
ammo-p.-

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

CROSSFLOW PERMEATION OF VISCOUS AND
VISCOELASTIC LIQUIDS THROUGH ARRAYS OF
CIRCULAR CYLINDERS

by

Craig A. Chmielewski

A DISSERTATION

Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Department of Chemical Engineering

1991

(you-- 3 /0,)§

ABSTRACT

CROSSFLOW PERNIEATION OF VISCOUS AND
VISCOELASTIC LIQUIDS THROUGH ARRAYS OF
CIRCULAR CYLINDERS

By

Craig A. Chmielewski

The ﬂow of viscous and viscoelastic liquids transverse to periodic arrays of
circular cylinders is studied. Two array geometries, square and hexagonal, are
examined; each having a void fraction of- 70 percent. Particular attention is devoted to
the inﬂuence of macromolecular conformation on the enhanced pressure drop across the
arrays. The ﬂow kinematics in both arrays is elucidated by laser'Doppler velocimetry
and streak photography. These techniques reveal a ﬂow transition from steady to
unsteady motion at Deborah numbers corresponding to the onset of enhanced pressure
drop. This result indieates that any attempt to predict the relative ﬂow resistance
increases observed with the viscoelastic ﬂuids must describe the transition to unsteady
ﬂow.

The elastic ﬂuids consist of four polyisobutylene (PiB) solutions; three non-shear
thinning and one shear thinning. The non-shear thinning liquids are 0.2 wt. 96 P13 in
polybutene solutions which differ only by the molecular weight of the solute. As these

solutions are 0—systems, the square of the degree of polymer extensibility associated with

Craig A. Chmielewski
each of these solutions is proportional to their molecular weight. This is consistent with
extensional viscosity measurements made by ﬁber spinning.

The initial departure from Darcy’s law for the non-shear thinning solutions is an
enhancement in ﬂow resistance and occurs at Deborah numbers of 0.80 and 0.35 for the
square and lwxagonal arrays respectively. These onset values are independent of
molecular weight. At a given Deborah number above the onset, the ﬂow resistance is
greater for higher molecular weights. Even after resealing to obtain the same onset
Deborah number for the two arrays, the relative ﬂow resistance is higher for the
hexagonal array than the square array. At large values of Deborah number the relative
ﬂow resistances in both arrays become independent of Deborah number. These
asymptotic values are shown to be proportional to the molecular weight and hence to the
square of the polymer extensibilities for these 0-systems.

The initial departure from Darcy’s law for the shear thinning solution is a
decrease in the relative ﬂow resistance followed by an increase. This increase occurs at

Deborah numbers of 4.2 and 3.0 for the square and hexagonal array respectively.

To my wife, Stephanie, with whose love and patience
has made this all possible.

iv

Acmomnnomms

I wish to extend special thanks to my thesis advisor, Dr. Krishnamurthy
Jayaraman, for his guidance and friendship throughout my graduate study. I also would
like to thank Dr. Charles A. Petty who always seemed available for an intellectual
discussion. I would like to thank the other members of my committee as well, for their
time and assistance: Dr. Lawrence T. Drnl, Dr. Dahsin Liu and Dr. David Yen.
Further thanks need to be extended the Department of Chemieal Engineering, the
Composite Materials and Structures Center and the State of Michigan without whose
ﬁnancial support this study could not have been undertaken.

I also need to acknowledge those who have provided special assistance to speciﬁc
portions of this study. These people include Drs. John Foss and Syed Ali for their
technieal and equipment support during the LDV studies, and Mike Rich and Brian Rook
of the CMSC laboratory for their support with the composite processing equipment.
Also, I would like to thank Terry Casey and Mike McLean of the D.E.R. machine shop
for their assistance in the construction of the permeation cells.

Finally, I would like to thank the many friends I have made over the last few
years who have made my stay at MSU so wonderful, especially Steve Reiken, Greg
Kudlac, Venkatesh Rae, Brent Larson, Nancy Losure, Kevin Nichols, Mike Bly, Shri

Iyer and Javad Kalantar.

TABLE OF CONTENTS

LIST OF TABLES ......................
LIST OF FIGURES ,. . . ....................
NOMENCLATURE ......................
Chapter
1. INTRODUCTION ................
1.1 Motivation ................
1.2 Newtonian liquid ﬂow in porous media . . . .
1.3 Non-Newtonian ﬂow in random media .......
1.4 Corrugated tube models .............
1.5 Planar elongation .............
1.6 Dissertation outline .............

2. THE EFFECT OF POLYMER EXTENSIBILITY
ON THE FLOW OF POLYMER SOLUTIONS
THROUGH CYLINDER ARRAYS .......

2.1 Summary ...................
2.2 Introduction ................
2.3 Experimental ................
2.3.1 Materials .............
2.3.2 Shear ﬂow properties ..........
2.4 Fiber spinning ................
2.4.1 Extensional ﬂow apparatus .......
2.4.2 Extensional ﬂow results and discussion .
2.5 Crossﬂow through cylinder arrays .......
2.5.1 Permeation apparatus ..........
2.5.2 Flow unsteadiness ..........

vi

H

:‘OQON'Fv—t

2.5.3 Polymer degradation ............
2.5.4 Permeation of Newtonian ﬂuids .......

2.5.5 Permeation of viscoelastic liquids .......
2.6 Conclusions . ..................

THE KINEMATICS OF VISCOUS AND
VISCOELASTIC LIQUID FLOWS WITHIN
ARRAYS OF CIRCULAR CYLINDERS .........

3.1 Summary .....................
3.2 Introduction . ..................

3.3 Laser Doppler velocimetry ............
3.3.1 Basic principles ............
3.3.2 Experimental LDV system .........
3.3.3 LDV measurement difﬁculties ......
3.4 Experimental ..................
3.4.1 Test ﬂuids ................
3.4.2 Apparatus . ...............
3.4.3 Streak photography ............
3.5 Validation of experimental technique .........
3.5.1 Stokes ﬂow simulation .........
3.5.2 Newtonian liquid ﬂow visualization ......
3.5.3 Newtonian liquid LDV measurements . . .
3.6 Elasitc liquid results and discussion .........
3.6.1 Flow visualization on the PIB/PB liquid
3.6.2 LDV measurements on the PIE/PB liquid . . .
3.6.3 LDV measurements on the PIB/decalin liquid . .
3.6.4 Discussion of ﬂow unsteadiness . ......
3.7 Conclusions . ..................

THE DEGRADATION OF POLYMER
SOLUTIONS FLOWING THROUGH
ARRAYS OF CIRCULAR CYLINDERS .........

4.1 Summary .....................
4.2 Introduction ...................

4.3 Experimental ..................
4.3.1 Materials ................
4.3.2 Rheological properties . ..........
4.3.3 Apparatus . ...............
4.3.4 Wall and end effects ............
4.4 Results and discussion ...............

Page

47
56

57

57
58
61
61
67

71
71
71
74
75
75
75
78
78
78
89

102

APPENDIX

4.4.1 Transverse permeability of Newtonian liquid . .
4.4.2 Departure from Darcy’s law for elastic liquids - .
4.4.3 Degradation of elastic liquids .........
4.4.4 Analysis of chain extension in arrays .....
4.4.5 Chain extension modeling .........
4.5 Conclusions ..................

VISCOSITY EFFECTS IN THE
PROCESSING OF COMPOSITE
PREPREG BY HOT MELT IMPREGNATION ......

5.1 Summary .....................
5.2 Introduction ..................
5.3 Experimental ..................
5.4 Resin and mass fraction and distribution . ......
5.5 Mathematieal modeling ...............

5.5.1 Tow consolidation models .........

5.5.2 Impregnation model ............

5.5.3 Solution procedure ............
5.6 Calculation results and discussion .........
5.7 Conclusions ..................

WEDGE DIFJFIBER IMPREGNATION
FORTRAN PROGRAM ...............

A.1 Flowchart . ..................

A.2 Subprogram summary ...............
A.3 Program .....................

viii

Page

115
117
121
125

136

LIST OF TABLES

Page
Table 2.1 Properties of the test ﬂuids at 20 °C. ........... 17
Table 2.2 Operating parameters in ﬁber spinning experiments ...... 29
Table 3.1 Components of the laser Doppler velocimetry system shown
in Figure 3.1. .................... 63
Table 4.1 Properties of the fresh test ﬂuids. ........... 109
Table 4.2 Properties of elastic liquids before and after crossﬂow
runs in different arrays .................. 123

ix

1.1

1.2

2.1 (a)
2.1 (b)
2.2 (a)

2.2 (b)
2.3 (a)
2.3 (b)
2.4

2.5 (a)

2.5 (b)

2.6

2.7

LIST OF FIGURES

The evolution of the viscosity and relaxation time

ofa5minuteepoxyresincuringatroomtemperature. . . .

The relative extensional viscosity versus De of a FENE

dumbbell liquid in a steady planar extensional ﬂow. . . .
Steady shear viscosity of the PIB in PB solutions. . . .

Steady shear viscosity of the PIB in deealin solution. . . .

Relative viscosities of the PIB in PB solutions as a
function of solution concentration.

Intrinsic viscosities of the PIB in PB solutions .......
Dynamic storage moduli of the PIB in PB solutions. . . .

Dynamic moduliofthePIBindecalinsolution. . . .

Schematic of ﬁber spinning apparatus. .........

Transient extensional viscosities as measured by
ﬁber spinning for the PIB in PB solutions.

Apparent Trouton ratios, accounting for the polymer
contribution, based on ﬁber spinning measurements

takenaté~15s“forthePIBinPBsolutions. . . .

Apparent Trouton ratios based on ﬁber spinning data of
the PIB in decalin solution where the vlues of II are

determinedaty=ﬁé ................

10
20

21

22
23
25
26

27

32

33

36

37

Rem
2.8 Schematic of the permeation cells (all dimensions in cm).

2.9 Geometry of the circular cylinders (radius = 0.238 cm)
(a) square pitch (b) hexagonal pitch. ............

2.10 (a) Pressure traces for the 2.11 x 10‘ molecular weight
PIB/PB solution in the hexagonal array at De = 2.6
(Q = 81 cm’ls, upstream P = 29.2 x 10‘ Pa, downstream
P = 1.1 x 10" Pa). ..................

2.10 (b) Pressure traces for the 2.11x10‘ molecular weight
PIB/PB solution in the square array at De= 2. 6
(Q = 113 cm’ls, upstream P = 17. 8 x 10‘ Pa, downstream
P = 1.7 x10‘ Pa). ..................

2.11 Friction factor vs. Reynolds number of two
Newtonian PB liquids compared to the theoretical
prediction of Sangani and Acrivos (1982) in square
and hexagonal arrays (closed symbols represent the
higher molecular weight PB). ...............

2.12 (a) Friction factor vs. Reynolds number of the PIB in
PB solutions ﬂowing through the square array ..........

2.12 (b) Friction factor vs. Reynolds number of the PIB in
PB solutions ﬂowing through the hexagonal array. ......

2.13 (a) Flow resistance vs. Deborah number for the PIB in
PB solutions ﬂowing through the square array ..........

2.13 (b) Flow resistance vs. Deborah number for the PIB in
PB solutions ﬂowing through the hexagonal array. ......

2.14 The high Deborah number ﬂow resistance asymptotes as
a function of molecular weight for both the square
and hexagonal pitch arrays. ...............

2.15 Flow resistance vs. Deborah number for the PIB in

decalin solution ﬂowing through both the square
and hexagonal array. ..................

xi

Page
39

48

50

3.1

3.2

3.3
3.4

3.5

3.6

3.7

3.8

3.9 (a) .

3.9 (b)

‘ array along y = 1.74. ...............

3.10

3.10 (0

Schematic of laser Doppler velocimeter system

(components are identiﬁed in Table 3.1). ......

Measuring volume and fringe pattern formed by
the beam intersection. The beams intersect in

the x-z plane and the z axis follows the bisector

of the angle of beam crossing (Dabir, 1983). ......

A typical signal from the photomultiplier (Dabir, 1983).

Steady shear viscosity of a 0.25 % PIB/PB solution

at25 °C. .....................

Linear viscoelastic properties of a 0.25 % PIB/PB

solution at 25 °C ...................

Geometric domains used for the Stokes ﬂow simulations

(a) square array and (b) hexagonal array. ......

Streamline output from the Stokes ﬂow simulations

(a) square array and (b) hexagonal array. ......

Streak photographs for a Newtonian polybutene liquid
in the cylinder arrays at a void fraction of 70 percent
(ﬂow from right to left) (a) square array, Re = 0.027

and (b) hexagonal array, Re = 0.013. .........

LDV measurements and the Stokes ﬂow prediction
for the inelastic Newtonian ﬂuid in the square

array along y = 1.62. ...............

LDV measurements and the Stokes ﬂow prediction
fortheinelasticNewtonianﬂuidinthehexagonal

Streak photographs for the viscoelastic liquid in the
square array having a void fraction of 70 percent
(ﬂow from right to left) (a) Re = 0.039, De = 0.16

(b) Re = 0.19, Dc = 0.80. ............

Re=0.25,De=1.09(d)Re=0.36,De=1.44. . . .

xii

73

Figure
3.10 (e)

3.11

3.12

3.13 (a)

3.13 (b)

3.14

3.15 (a)

3.15 (b)

3.16

4.1

4.2

4.3

Re=0.48,Dc=1.9l(0Re=0.48,De=1.91. . . . .

Streak photographs for the viscoelastic liquid in the
hexagonal array having a void fraction of 70 percent
(ﬂow from light to left) (a) Re = 0.015, De = 0.059

(b) Re = 0.069, De = 0.28. .............

3.11 (c)Re = 0.11, De = 0.43 ((1) Re = 0.14, De = 0.56.

Relative ﬂow resistance of a 0.25 % PIB/PB elastic liquid

in both square and hexagonal cylinder arrays. .......

LDV measurements and the Stokes ﬂow prediction for the
0.25 96 PIB/PB elastic liquid in the square array along

y = 1.62. ......................

IDV measurements and the Stokes ﬂow prediction for the
0.25 % PIB/PB elastic liquid in the hexagonal array

along y = 1.74. ...................

Probability distribution function of the elastic ﬂuid’s
velocities measured in the hexagonal array at x = -1.58

andy = 1.74 (a)De = 0.085 mm = 0.36 ........

LDV measurements and the Stems ﬂow prediction for the
PIB/Decalin liquid in the square array along y = 1.62.

LDV measurements and the Stokes ﬂow prediction for the
PIB/Deealin liquid in the hexagonal array along y = 1.74.

Growth of the ﬂow instability in the hexagonal array for
the 0.25 % PIBIPB elastic liquid as measured by LDV at

x = -1.58 and y = 1.74.. ................

Steady shear viscosity ofthe fresh test ﬂuids. .......

Dynamicmoduliandsteadyshearviscosityofthe

fresh Ml ﬂuid. ...................

(A) Schematic of experimental apparatus and
(B) a cross-sectional view of the permeability

cell(alldimensionsareincentimeters) ...........

xiii

Page
85

87

88

91

93

94

.95

..100
..107

. .110

. .112

4.4

4.5 ,

4.6

4.7

4.8

4.9

4.10

4.11

4.12

4.13

5.1

5.2

Page
Cylinder array geometries: (A) rectangular pitch array;
(B) triangular pitch array (radius = 0.159 cm and
porosity = 0.704). .................. 114
Mean pressure gradient vs. superﬁcial velocity for
the Newtonian ﬂuid. .................. 116
Friction factor vs. Reynolds number for the M1 ﬂuid
ﬂowing through the triangular pitch array. ......... 119
Friction factor vs. Reynolds number for the M1 ﬂuid
ﬂowing through the rectangular pitch array. ......... 120
Flow resistance of the M1 ﬂuid relative to the
Newtonian value vs. Deborah number. ............ 122
Streamlines computed by ﬁnite element simulation
of Stems ﬂow through (A) the rectangular pitch
array and (B) the triangular pitch array (i and
ii denote selected streamline segments) ............. 126
Velocity proﬁles along selected streamline segments
(see Figure 4.9). ..................... 128
Strainrateproﬁlesalongselectedstreamlinesegments
(see Figure 4.9). ..................... 129
Computed stretch ratio of a polymer chain along selected
streamline segments: De = 0.4 ................ 131
Computed stretch ratio of a polymer chain along selected
streamline segments: De = 1.0 ................ 132
Hot melt impregnation (a) schematic of the hot melt
prepregger (b) schematic of the resin pot. ......... 138
Resinmassfractiondependenceontowvelocityand
resin viscosity. ..................... 140

xiv

5.3

5.4

5.5

5.6

5.7

5.7

5.7

5.8

5.8 '

5.9

A.l

Two 12k AS4 earbon ﬁber tows impregnated with a
DG-A/mPDA epoxy resin system under the following
processing conditions: tow velocity = 18.3 calls

(a) resin viscosity = 4.5 Pa-s (b) resin viscosity

= 0.1Pa-s. ....................

An illustration of the straightening of a portion of

a bent ﬁber isolated from a network of ﬁbers;

(a) the buckled portion of ﬁber under consideration

in an unstressed ﬁber network; 0)) as the ﬁber network
is placed under tension, the ﬁber portion of interest
takes up some of the load; (c) the result is that the ﬁber

straightens ...... A .................

Model non-linear tow consolidation behavior as a result of

tensionontheﬁbertow. ................

Wedge die geometry. . ................

(a) Calculated pressure proﬁles for several different

values of the effective tensile modulus. ..........

(b) Calculated impregnation proﬁles for several different

values of the effective tensile modulus. ..........

(c) Calculated ﬂow rate proﬁles for several different

values of the effective tensile modulus. ..........

(a)Calcu1atedresindistributioninamodelﬁbert0w. . . . .

(b)Calcu1atedresinmassfractionofamodelﬁbertow.
Relationship between the partially impregnated resin tow

tensile modulus and the resin viscosity. ..........

Schematic of the computer algorithm. ..........

. .147

..161

163

..164

165

.. 167

..174

cr an

O

agrees

p

"’5“

KR)

NOMENCLATURE

cylinder radius

Mark-Houwink exponent

beam width

concentration

fringe spacing

rate of deformation tensor
diameter

Deborah number

Deborah number deﬁned with the packing number
draw ratio

beam modulus

friction factor, foeal length

force law for the FENE dumbbell
applied force on beam

ﬁber tension

superﬁcial mass flux

storage modulus

xvi

G"

K11 9K22 rK33

KI

loss modulus

beam width

die shape

coating thickness

Kozeny constant

permeability tensor

diagonal components of the permeability tensor

transverse permeability

'Mark-Houwink pre-exponential factor

cylinder array length

beam length

extensibility parameter, die length
liquid ﬁber length

viscosity averaged molecular weight
weight averaged molecular weight
cycles per burst

Capillary number

atmospheric pressure

pressure

compressive load

ﬂow rate

ﬁber radius

r,R

VO’VO

end-to-end distance of a polymer molecule
Reynolds number

impregnation front

ﬁber tow proﬁle

effective impregnation

time 1

temperature

rr-component of the total stress tensor
22-00th of the total stress tensor
normal velocity component into the ﬁber tow
ﬁber tow velocity

velocity

superﬁcial velocity

x-component of the velocity

y-component of the velocity

available ﬁber volume fraction

ﬁber volume fraction

ﬁber volume fraction based on the entire die cross-sectional area
initial ﬁber volume fraction

die width

die sap

x coordinate direction

x .0.

’Io
’ln
’ln
’11:
ll:
’7-

[n].

y coordinate direction

2 coordinate direction
W

packing factor

shear rate

resin-ﬁber surface tension

capillary pressure

void fraction

extension rate

beam strain

shear viscomty

zero shear viscosity

extensional viscosity

transient extensional viscosity

average extensional viscosity

relative viscosity

solvent shear viscosity

intrinsic viscosity

die half angle

half angle of laser beam crossing

spring constant (compressive loads)

spring constant (tensile loads)

xix

kt

”n

laser light wavelength
relaxation time
Doppler frequency
shift frequency
nasty

contact angle
frequency

Sperm
signiﬁes averaging of "

superscript asterisk signiﬁes x is a dimensional quantity

subscript zero signiﬁes x is an inlet value

Chapter 1

INTRODUCTION

1.1 Matizattim

The ﬂow of rheologically complex ﬂuids in random porous media has long been
of interest to the oil industry for use in oil recovery processes. With conventional
recovery methods, a large amount of oil often remains in the ground after regular
production ceases, requiring other innovative techniques to extract the residual oil.
These ancillary techniques fall under the heading of ”tertiary oil recovery”. Tertiary oil
recovery processes typically involve pumping dilute polymer solutions in wells to
displace inaccessible oil from porous rock. Small amounts of high molecular weight
polymers are added to increase the viscosity of the displacing ﬂuid and prevent viscous
ﬁngering during permeation. However, the ﬂuid rheology becomes complex,
complicating the recovery process at high permeation rates.

Themateﬁalprocessingmdumyalsohasgareratedanecdmunderstandtheﬂow
ofviscousandviscoelasticliquidsthroughporousmedia, suchasﬁbermats, inorderto
produce lower cost and higher quality composite materials. Traditionally, the challenge
had been to wet out bundles of ten to twenty micron ﬁbers with viscous epoxy resins.

1

2
More recently this problem has been expanded to wetting out ﬁber bundles with

viscoelastic ﬂuids. The complex rheology of high molecular weight materials and fast
curing epoxy resins complicate composite manufacturing processes by introducing elastic
effects. The magnitude of these effects may vary, and can be quantiﬁed by a
dimensionless group named the Deborah number, De; the ratio of the characteristic ﬂuid
time scale to the ﬂow or processing time scale. For example, impregnation and
consolidation ﬂows of rubber modiﬁed resins and thermoplastics, having time constants
in the range of 0.1 to 10 seconds, achieve Deborah numbers of order one as they ﬂow
past 10 micron ﬁbers at velocities as slow as 0.0001 to 0.01 cm/s.

Continued growth in the use of rapidly curing thermosetting systems for liquid
molding operations also provides impetus for studying viscoelastic ﬂuid effects in
composite processes. Figure 1.1 shows the time evolution of the viscosity and relaxation
time of a 5 minute epoxy, undergoing a thermosetting reaction at room temperature. The
initial magnitude of the relaxation time is relatively low, but quickly develops so that
over the last 50 percent of the epoxy’s ﬂuid life the liquid is considerably elastic. This
sort of rheological complexity affects the time and pressures needed to ensure the
complete saturation of a ﬁber preform.

Thus, complex ﬂuid rheology must be accounted for in order to obtain a basic
understanding of the ﬂow of high molecular weight (or growing molecular weight)
materials in porous media. This is particularly a concern to the composite industry,

where the ability to move these advanced materials into high volume commercial use

 

 

 

 

* .1000
154-04 . / E
Tam . - 27 'C _
P / _ m
r- 9..
o
,_, 1000 g
9 a:
g 8
H 71.0 _1
100 : _
a? . §
m h
8 ' a“
9 ' H
> 10
U " ‘I
1 ' 4 e
0 Y

Time [min]

Figural 1 'Iheevolutionoftheviscosityandrelaxationtimeofasminuteepoxyresin
reectingatroom temperature.

4
depends in part on reducing processing time. A fundamental approach to the solution of

this problem entails the study of the dynamics of liquids, characterized rheologically
under both shear and extension, ﬂowing through uniform arrays of circular cylinders.

11W

The earliest scientiﬁc work on low Reynolds number ﬂows of Newtonian liquids
through isotropic porous media was done by H.P.G. Darcy (1856). He found that the
superﬁcial velocity of ﬂuid through a porous medium was proportional to the pressure
gradient across the medium. This result is now known as Darcy’s law, and in its general

form is:
v, - .’T‘<Vp> . (1.1)

Here v.is the superﬁcial velocity, u is the shearviscosity, Kis thepermeability tensor,
deﬁned solely by the geometry of the media and <VP> is the mean pressure gradient
in the liquid. For random media, such as packed beds of spheres, the permeability
tensor is isotropic. However, for the uniform arrays of circular cylinders studied here,
four components of the tensor in general are necessary to describe the permeability: K“,
K,,, Km and K3,. In addition, for cylinder arrays with rotational symmetry, such as

square and hexagonal pitches, there are only two unknown components, the transverse

and axial permeabilities.

The study of the ﬂow resistance offered by arrays of cylinders has its origins in
the design of tubular heat exchangers. The need for more efﬁcient heat transfer
equipment continues today, and has prompted the theoretical study of low Reynolds
number viscous ﬂows through periodic arrays of cylinders by Sangani and Acrivos
(1982). They calculate the drag of a Newtonian liquid on cylinders arranged in square
and hexagonal pitches over a wide range of loadings. For void fractions greater than 70
percent, both array geometries offer the same resistance to ﬂow; with void fractions of
50 percent or lower, the square array offers more resistance to the ﬂow of a Newtonian
liquid.

Recently, the composite industry has motivated further study of permeation
through ﬁber preforrns. Adams et al. (1986) have developed a planar ﬂow technique to
studythepermeabilitiesofwoven ﬁber mats. Intheirexperimentsanepoxyresinwas
injected under a constant pressure into a variety of ﬁber mats, each having various weave
patterns and ranging in void fractions from so to 91 percent. The epoxy was forced to
ﬂowradiallyintheplaneofapreform, andtheprogress oftheﬂowfrontwas monitored
over time. For the case of isotropic mats, the ﬂow front was circular, indicating a
radially uniform in-plane permeability which could be backed out by an unsteady material
balance and Darcy’s law. In the ease ofanisotropic mats, however, the ﬂow front was
elliptic, but also could be predicted using an unsteady material balance and Darcy’s law.
The in-plane permeability was not constant, however, and had to be described by a linear

combination of the permeabilities measured along the principal axes of ﬂow (the

6
directions of maximum and minimum flow). More recently, Adams and Rebenfeld

(l991a,b) have used this same technique to study the in-plane flow through multilayer
fabric assemblies, and found that the overall permeability differs from that of the
individual constituent layers. Furthermore, they found that flow transverse to the plane

was important in maintaining a macroscopically uniform flow front.

1.3 W

One of the earliest studies of non-Newtonian flow through porous media was that
of Sadowski and Bird (l965a,b). They studied the flow of shear thinning aqueous
polymer solutions through random media. These solutions were virtually inelastic; the
maximum time constant being 0.07 seconds. They correlated their data by modifying
Darcy’s law to include the fluid’s shear thinning viscous behavior. This was done by
replacing the constant viscosity in Darcy’s law with an appropriate generalized
Newtonian ﬂuid model-inthiscasetheEllismodel-whoseparameters werechosenby
ﬁttingtheshearviscositydata. Thismethod described thepressuredrop-ﬂowratedata
well, and has been used in subsequent analyses of the permeation of weakly elastic, shear
thinning liquids through random porous media ((12 Christopher and Middleman, 1965).
It was not until the work of Marshall and Metzner (1967) that sufﬁciently high nominal
bedstrainrateswerereached, resultinginﬂowresistancesuptotentimestheNewtonian
value. Theseresearchersargued thattheenhancedpressuredropwastheconsequence

of molecular extension due to the converging-diverging nature of the packed bed; a

7
featurethecapillarymodelsdonotcontain.

1-4 W

The inadequacy of the combination Darcy law/ generalized Newtonian ﬂuid models
to predict the large increases in ﬂow resistances of polymer solutions ﬂowing through
randommediahasleadtonewefforts-toformulatemodelswhichtakeintoaccountthe
converging-diverging nature of porous media. The model which has garnered much
attention lately is the corrugated tube which is an axisymmetric tube with a sinusoidally
varying radius. Over the past twenty years several experimental and numerical
investigations concerning the ﬂow of elastic liquids in corrugated tubes has been
undertalmn (cf Dodson er al., l97l).

James er al. (1990) report experimental measurements ofpressure drop in slow
ﬂow of an elastic liquid through a corrugated tube. Their results show little deviation
fromNewtonianbehavior,evenatDe = 3,andagreeverywellwiththenumerical
predictions of Pilitsis and Beris‘ (1989) and Burdette er al. (1989). The more recent
numerical work includes calculations by Pilitsis and Beris (1989), who use an upper
convected Maxwell model and a pseudospectral/ﬁnite difference method, and Burdette
er al. (1989), who also use an upper convected Maxde model but an explicitly elliptic
momentum equation formulation. The results of these two numerical experiments are
inexcellentagreement, evenatDe =10, and, liketheexperiments, predictvirtuallyno

variationintheﬂowresistancefromtheNewtonianvalue. Onlythroughtheinclusion

8
of inertia have Pilitsis and Beris (1991) been able to predict a substantial enhancement

in the ﬂow resistance. These results, both experimental and numerical, indicate that the
conugatedtubemodelisnotadequatefordescribingtheﬂowofviscoelasﬁc ﬂuidsin
porous media.

A deﬁciency of the corrugated tube model is its lack of stagnation points (of
course, near the tube wall the ﬂuid moves slowly as a result of the no slip condition).
The term stagnation point will be used here not only to refer to a point of zero velocity,
but also to a point situated in a ﬂow region containing large extensional deformation
rates. A macromolecule, ﬂowing into the neighborhood of a stagnation point, will
experience high residence times while simultaneously being stretched by the extensional
ﬂow. Hence, a stagnation point within the ﬂow ensures the existence of the necessary
criteria for a coiled macromolecule to undergo elongation: large residence times and
large extension rates.

Stagnation points abound in porous media. For example, each sphere comprising
the model random media of the above experiments contains two stagnation points, each
located at the poles. Neglecting the complications resulting from sphere-sphere contacts,
theﬂuiddirectlyupstreamofasphereexperiencesabiaxialextensionwhiletheﬂuid
directly downstream of a sphere undergoes uniaxial extension. A similar situation occurs
intwodimensionsfortheﬂowothuidsu'ansversetocircularcylinders. Acylinder
alsohastwostagnationpoints;againeachlocatedatthepoles. Becausctheﬂowistwo
dimmsional,meﬂuidbomupsueunmddownmeamofmesmgnaﬁmpommexperiurces

planar elongation.

1.5 Wan

To demonstrate the effect an extensional ﬂow has on ﬂow dynamics, we shall
examine the stress developed in a model viscoelastic liquid undergoing a pure
elongational ﬂow; the relevant kinematics for this study being planar elongation. The

velocity components for a steady planar elongational ﬂow are given by,

vJr = xé (1-2)

V, = -yé , (1.3)

whereéisﬂreextensionmteandxandyrepresentperpendicﬂarCanesiancoordinate

directions. The rate of deformation tensor for this ﬂow is represented by the matrix, D,

,, . .[g 31} . (1...,

For this example, a dilute solution of finitely extensible non-linear elastic (FENE)
durirbbells will be used as a ”realistic" representation of a dilute polymer solution. In
thepresentcontext,theFENEdumbbellmodelisrealisticinthesensethatthetotal
extension of its macromolecular components, the dumbbells, is limited. This prevents

the stress from becoming inﬁnitely large in regions of high extension rates. Figure 1.2

10

 

 

 

 

>\

’5, 1000.0
0

0

92

>

'6 100.0
C

O

'63

C

3 10.0
X

U

0

..>. 1 0
E o
0

(K

0.1 I V 1 1 V V V VI f ﬁ V V U U V 1
0.1 1.0 10.0

Deborah Number

1113111212 TherelativeextensionalviscosityversusDeofaFENBdumbbenliquidin
asteadyplanarextensionalﬂow.

l 1
shows the steady planar extensional viscosity behavior with Deborah number of Chilcott

and Rallison’s (1988) FENE dumbbell liquid. The important features of this model
prediction are that over a very small range of Deborah numbers the extensional viscosity
increases appreciably over the Newtonian value, and that for sufﬁciently high Deborah
numbers the extensional viscosity becomes independent of the Deborah number. Based
on this prediction, it is reasonable to expect that at an onset Deborah number the relative
ﬂow resistance of viscoelastic liquids, ﬂowing through arrays of cylinders (and for that
matter, through any porous media peppered with stagnation points), will increase and

then eventually become independent of the Deborah number.

1.6 W

This study is a fundamental investigation of the ﬂow of viscous and viscoelastic
liquids through arrays of circular cylinders. Chapter two focuses on the problem of how
geometry and ﬂuid rheology affect the ﬂow resistance through model arrays of circular
cylinders. Rheologically characterized ﬂuids in both shear and extension are used to
elicit what ﬂuid and geometric parameters are important in predicting the dynamic
behavior of ﬂuids in cylinder arrays. In chapter three laser Doppler velocimetry and
streak photography are used to examine how the kinematics of the cylinder array ﬂows
are affected by the ﬂuid rheology. Chapter four presents experimental data showing the
effects of array geometry on macromolecular chain scission of a high molecular weight

polymer solution ﬂowing through rectangular array and in a hexagonal arrays. Chain

12
extension calculations based on the Stokes ﬂow ﬁeld are presented to support the

experimental observations. Finally, whereas the previous chapters consider the effect of
ﬂuid elasticity on the permeation of liquids through cylinder arrays, chapter ﬁve
considers the effects of cylinder array compliance on crossﬂow permeation. This is done
through an analysis of the hot melt impregnation process which is used in the
manufacture of composite prepreg. In this chapter it is shown that ﬁber bundle elasticity

has a profound effect on the permeation of viscous liquids into ﬁber arrays.

Chapter 2

THE EFFECT OF POLYMER EXTENSIBILITY ON
THE FLOW OF POLYMER SOLUTIONS THROUGH

CYLINDER ARRAYS

2.1 Smut

The effects of ﬂuid rheology on low Reynolds number ﬂows transverse to
periodic arrays of circular cylinders have been investigated with several solutions of
polyisobutylene. Care was taken to avoid degradation of the polymer during the
measurements. These solutions were theologically characterized in both shear and
extension. Three of the solutions are dilute solutions of different molecular weight
polyisobutylenes in polybutene at the same concentration. These are o-systems at room
temperature and have a constant shear viscosity over strain rates up to 10 3". Fiber
spinning of these solutions indicates that the apparent Trouton ratio at an average stretch
rate of 15 s" is proportional to the molecular weight. This is consistent with predictions
of FENE dumbbell models for higher elongation rates (cf Chilcott and Rallison, 1988;

Biller et al. , 1986). A fourth solution of polyisobutylene in decalin was used to evaluate

13

14
the ﬂow resistance for shear thinning solutions.

The resistance to ﬂow of the non-shear thinning solutions in both square and
hexagonal pitch arrays is above the Newtonian value at onset Deborah numbers of 0.80
and 0.35 respectively. These onset values are independent of solute molecular weight.
Higher molecular weight ﬂuids produce higher ﬂow resistances relative to the Newtonian
value for Deborah numbers greater than the onset value. At large Deborah numbers (De
> > 1) the relative ﬂow resistances in both arrays become independent of Deborah
number, and scale linearly with the molecular weight. The asymptotic value of the
resistance ratio is consistently higher for the hexagonal array than for the square array
for the same molecular weight. This is shown to be a result of these transverse ﬂows

being dominated by planar extension at high Deborah numbers.

mm

It is well recognized that the dynamics of shear thinning polymer solutions may
be affected by the ﬁnite extensibility of its macromolecular components (cf. , Christiansen
and Bird, 1977/1978). In spite of this, the effect of varying polymer extensibility on the
ﬂuid dynamics of non-shear thinning elastic liquids, the so-called Boger liquids (Boger,
1977/ 1978), has not adequately been explored. A reason for this may be that the
polybutene based Boger liquids have been shown to behave in both shear and extension
as a dilute solution of inﬁnitely extensible, linear dumbbells at low to moderate

deformation rates (see Prilutski er al., 1983 and Sridhar er al., 1986 respectively). In

15
strong ﬂows, however, where the Deborah number (De) is much greater than one, this

model breaks down because a linear dumbbell will extend without bound in such a ﬂow
(Rallison and I-Iinch, 1988). This may be one reason that numerical simulations using
Oldroyd type constitutive models fail quantitatively, and sometimes even qualitatively,
to account for important ﬂow phenomena observed for Boger ﬂuids undergoing complex
ﬂows. Only recently have researchers begun to explore other models in order to better
understand the physics behind the ﬂow properties of polymer solutions. For example,
Chai and Yeow (1990) use a multiple relaxation time constitutive model (KBKZ model)
to describe the ﬂow of a Boger ﬂuid in a gravity drawn jet. They found better
agreement with experimental data using the KBKZ model than with the Oldroyd-B
model.

In the present investigation the effect of varying ﬁnite extensibility of the polymer
on the dynamics of complex ﬂows of Boger liquids will be examined. This study is
partially motivated by the work of Chmielewski er al. (1990a) who observed opposite
trends with Deborah number in the relative drag on spheres, translating in corn syrup
based polyacrylarnide solutions versus polybutene based polyisobutylene solutions. This
difference was attributed to differences in the extensions of the polymers from
equilibrium. Polyacrylamide molecules tend to be relatively elongated in aqueous
solutions at equilibrium while polyisobutylene molecules in polybutene tend to be
relatively coiled at equilibrium.

In a previous investigation, Chmielewski er al. (1990b) reported that the pressure

drop in ﬂow transverse to arrays of circular cylinders is much greater for an elastic

16
Boger liquid than the pressure drop for Newtonian liquids of equivalent viscosity; similar

to that of polymer solutions in isotropic porous media. The focus of their study,
however, was on how differences in the extensional ﬂow ﬁeld between hexagonal and
rectangular pitch geometries affected polymer extension, and in turn molecular chain
scission.

In this work three Boger ﬂuids are prepared with varying degrees of polymer
extensibility in order to understand its effect on the ﬂow resistance of these ﬂuids in
cylinder arrays. Two array geometries are studied, square and hexagonal, both having
a void fraction of 70 percent. Unsteady extensional viscosity data, obtained by ﬁber
spinning, is correlated with high Deborah number ﬂow resistance asymptotes found in
both array geometries via molecular ﬁnite extensibility. A fourth, more concentrated
polymer solution is also examined in order to compare the effects of shear thinning on

ﬁber spinning and cross ﬂow resistance.

23W

2.3.1 Materials

Three classes of liquids were used in this study: two Newtonian liquids of widely
different viscosities, three non-shear thinning, elastic liquids which differed only by their
solute molecular weights, and a highly shear thinning, viscoelastic polymer solution. A
comparison of some of the material properties of these ﬂuids can be found in Table 2.1.

The two Newtonian liquids were 610 and 1290 weight average molecular weight

17

Table 2.1 Properties of the test ﬂuids at 20 °C.

FLUID

   

NEWTONIAN

 

 

 

 

   

Polybutene, H25

 

      

 

Polybutene, H300
7.45 % Kerosene

VISCOELASTIC
Boger Liquids

Solvent

 

 

 

 

 

93.0 96 PB, H25
7.0 96 Kerosene

 

Sohrtiom

 

0.20 % PIB L-80

0.90x10‘ "

 

0.20 96 PIB L-lOO

1.25x10‘ "

 

 

0.20 96 PIBL-14O

VISCOELASTIC
Shear Thinning

 

Solvent
Decalin

Solution
2 96 PIB B200

 

2.11x10‘ "

138.25

4.30x10‘

* Viscosity averaged molecular weight

** Zero shear viscosity

 

 

 

 

0.258

I 0.0027 I

 

18
polybutenes (PB) supplied by Amoco Chemical Company. A small amount of kerosene

(7.45 96) had to be added to the higher molecular weight grade in order to reduce its
viscosity to a manageable level. The three Boger solutions were all prepared by
dissolving an appropriate amount of polyisobutylene (PIB), supplied by Exxon Chemical
Company, into kerosene and then mixing the solution into Amoco’s 610 molecular weight
(MW) grade PB for a ﬁnal composition by weight of 0.20 96 PIB, 7.00 % kerosene and
92.80 % PB. Three PIBs, having different molecular weights, were used: (1) Exxon’s
Vistanex L-80 with a viscosity. average molecular weight of 0.90 (1: 0.15) x 10‘, (2)
Vistanex L-100 with a viscosity average molecular weight of 1.25 (:l: 0.19) x 10‘ and
(3) Vistanex L-l40 with a viscosity average molecular weight of 2.11 (:1: 0.24) x 10‘.
The shear thinning, viscoelastic liquid was a 2 96 solution of PIB (BASF B200)
in a cis and trans decalin mixture. This ﬂuid was supplied by Professor Walters, and its
preparation and composition were identical to the D1 liquid used in the Second Normal

Stress Diﬁ'erence Projea (Walters, 1983).

2.3.2 Shear ﬂow properties

Rheological measurements in steady and oscillatory shear were made on a
Rheometrics RFS-8400 Fluids Spectrometer at rates ranging between 0.1 and 100 s".
Thesetestswereperformedusinga0.02radianconeanda5cmdiameterplate.
Measurements made on the Newtonian and Boger liquids were performed at 10, 20 and
30 °C, while tests made on the remaining liquid were performed at 0, 20 and 40 °C.

Attempts to study the PIB/decalin solution at temperatures higher than 50 °C were

l9
unsuccessful due to enhanced solvent evaporation. Figure 2.1 (a) shows the solvent

viscosity, 1),, and the steady shear viscosity, 1;, of the three PBIPIB solutions as a
function of the shear rate. None of these liquids were signiﬁcantly shear thinning over
the range of shear rates tested. On the other hand, the shear viscosity of the PIB/decalin
solution was constant only up to 0.1 s", and then decreased with increasing shear rates
(see Fla-Ire 2.1 (b))-

Intrinsic viscosity, [1)],, measurements were made on the PBIPIB solutions to
assess the thermbdynamic solvent quality near 20 °C. This was accomplished by
calculating the relative viscosities, 11,, of three different concentrations of each of the
PBIPIB solutions. A linear extrapolation of ln(11,,)/c to inﬁnite dilution was made to
determine [1;], (see Figure 2.2 (a)). In Figure 2.2 (b) the relationship between [11], and

M, is shown to follow the Mark-Houwink relation,

[n], = K’M,‘ . (2.1)

The Mark-Houwink exponent, a, is approximately 0.5, indicating that the PBIPIB
solutions of this study were examined under theta conditions. The pre-exponential factor,
K’, is 0.40 cm’lg.

Figure 2.3 (a) shows the storage modulus, G', of the three PBIPIB solutions. A

relaxation time for each of the solutions (see Table 2.1) is calculated by,

20

 

 

 

 

10
. o 11,-0.90-10- 'r-zo'c
‘ e M'-1.25'10'
« n 11,-2.11-10-
4 I Solvent
7," 4
e 0 ' . O I g . .
O . . ' U I I I I
O. ' '
H d
0 O . . . .
p 0 O . . . . . C 0
' 0
II I I I I I I I I I I I t!
1 .T....... . ......., T.
0.1 1.0 10.0 100.0

Figure 2.1 (a) Steady shear viscosity of the PIB in PB solutions.

11(7) [Pa-S]

21

 

 

 

 

100.0—
0 T - 20 'C
0000000000000
°o
10.0-3 °o.°°
2 °o°
o
0.00
1.0::
1 °°o#
.1
0.1 1 . n"... . . ..m, .. ..n.
0 01 0.10 1.00 10.00 100.00
'5' [8"]

Figure 2.1 (b) Steady shear viscodty of the PIB in decalin solution.

22

 

T=20"C

|n(n.)/c [W 9]

1.04 o M «10.90ll 10‘
I M -1.25'10'
‘ n M =2.11'10‘

r I I
0.00 0.05 0.10 0.15 0.20
cone. [g/dl]

P
O
1:

 

 

 

FigureZJ (a) RelativeviscodtiesofthePIBinPB solutionsasafunctionofsolution
concentration.

[n] dI/g

10.0

dan‘

 

T=20°C

 

 

E+05 ' ' ”1364-06
Mel. Wt. [g/mol]

Figure 2.2 (b) Intrinsic viscosities ofthe PIB in PB solutions.

 

' 3:05.207

24

I
A . G . 2.2

 

The value of the storage modulus is taken from the low frequency region where G’ is
quadratic in frequency. Figure 2.3 (b) shows the dynamic moduli, at 20 °C, of the
PIB/decalin solution. The storage modulus of this ﬂuid also shows quadratic behavior

with frequency, but at much lower rates than the liquids of Figure 2.3 (a).

2-4 W

2.4.1 Extensional flow apparatus

The unsteady elongational ﬂow properties of the elastic ﬂuids were measured by
ﬁber spinning. The apparatus, shown schematically in Figure 2.4, is similar to that of
other researchers (see Hudson at al. , 1974, for example), and utilizes a bending beam
load cell to measure the force exerted by test ﬂuids on the capillary. Nitrogen pressure
is used to pump ﬂuid from a reservoir and through a 20 cm long stainless steel capillary
(0.238 cm I.D.). The capillary pivots at its upstream end via a miniature ball bearing
and is attached to the load cell at its downstream end.

The bending beam load cell consists of a 7.0 cm by 1.3 cm by 0.013 cm strip of
spring steel, actingasacantileverbeam. Thefreeend ofthebeamis attached totheend
of the capillary. At the ﬁxed end of the beam, four encapsulated foil strain gages

(Omega Engineering, model DYl l) are arranged to form a Wheatstone bridge. The

6' [Po]

 

 

100.000
o M_-0.90-10' T-20 'C 11
o M_-1.25-10- _;n
D -2.11'10' I.
10.000 M. g....
I':'
.e'o:
e'.o::.
1.000 .I'.::'
0.100 .' 0:.
0.010 ' '
0.001 r r r . r .1 r . . ”.1 r .n.
0. 10 10.0 100.0

Figure 2.3 (a) Dynamic storage moduli of the PIB in PB solutions.

 

0', G" [Pa]

26

 

100.00

° ‘3" T-20‘C

O G' ...,

0..
I...::OO....‘
10.00 .....:I II
..O::
O.....
1.00
O. ..

0.10 .'

.0

I
000 ‘ﬁ‘ T r U

 

 

Figure 2.3 (b) Dynamic moduli of the PIB in decalin solution.

27

 

 

 

 

 

 

 

 

 

 

 

 

 

maﬁa“ ‘
Um“
mac-awed -] m
J
*1
00"!
mm amt-1m
Tire-spawn
Qt
[— ‘1

 

Figure 2.4 Schematic of the ﬁber spinning apparatus.

28
bridge is powered with 10 volts and emits a millivolt signal proportional to the bending

strain in the beam. The beam strain, 15,, is proportional to both the applied force, F,,

and the moment arm, 1,, according to,

6 FBI,

5 a
" 2.12113

. (2-3)

where b is the width of the beam, h, is the beam thickness and E, is the beam modulus.
The load cell is sensitive to 0.1 mN and the signal output is linear over the range of
forces measured in this study (4 to 50 mN). The utility of this force measurement
method is not only its cost effectiveness, but it also provides a convenient way of
adjusting the sensitivity of the load cell to meet the needs of a particular application.
Steady extension of a liquid ﬁlament extruding from the capillary is maintained
by an adjustable speed take-up drum at the downstream end of the ﬁber. Still photograph
enlargements of an elongated liquid ﬁber are used to obtain ﬁlament diameter proﬁles.
Table 2.2 lists the operating parameters of the elongational ﬂow experiments performed

on each of the PIB test solutions.

2.4.2 Extensional flow results and discussion
The ﬁber spinning experiments were conducted under isothermal conditions and
such that gravity, surface tension and ﬂuid inertia were all very small compared to the

tensile force. The cross sectional average of the tensile stress at any axial location, 2,

29

Table 2.2 Operating parameters in ﬁber spinning experiments.

 

 

 

 

Boger unids

 

0.20 96 PIB L—80
0.20 96 PIB L-lOO

 

 

 

 

 

 

0.20 96 PIB L-140

Shear'l‘hinning

 

2%PIBB2001.

 

 

 

 

 

 

 

 

 

down the ﬁber was calculated by,

45’,
nD(z)3

 

T“(z) -r,,(z) = (2.4)

where T, and T, are components of the total stress tensor, F, is the tensile force and

D(z) is the diameter of the ﬁlament at location 2. The extensional strain rate will be

represented by the average, <é> ,

v(L,,) - v(0)

<e> =
LPS

 

. (2.5)

where LP, is the length of the ﬁber, and v(0) and van) are the axial velocities at the
capillary outlet and at the take-up drum respectively. Since the axial velocity proﬁles in
these spinning experiments are nearly linear, the local extensional strain rates at all axial
locations are approximately equal to the average extensional strain rate. Combining Eqs.

2.4 and 2.5, the transient extensional viscosity can be represented as,

T.‘(Z) - Tzr(z)

(t) (2.6)

 

n,‘(z) -

The ﬁber spinning results of the PBIPIB solutions are shown in Figures 2.5 (a)

31
and (b). In Figure 2.5 (a) the transient extensional viscosity is scaled with 311 and the

time is scaled with x,. This ﬁgure indicates a molecular weight effect on the transient
extensional viscosity independent of the dominant molecular time scale. An average
extensional viscosity between two points on the spin line may be calculated as suggested

by Mackay and Petrie (1989),

_ 4 15L”

'1: g -
zD(O)’v(0)1n(l‘-{1‘l) (2'7)
v(0)

 

The ratio of the extensional viscosity to the shear viscosity, accounting for the polymer
contribution, is plotted against molecular weight in Figure 2.5 (b) . Since these ﬂuids
had similar pre-shear histories and were spun at approximately the same extension rate
(5' ~ 15 s"), these results demonstrate that the apparent Trouton ratio scales linearly with
the molecular weight. This result can be anticipated by treating the Boger ﬂuids as a
dilute solution of ﬁnitely extensible non-linear elastic (FENE) dumbbells. At sufﬁciently
large Deborah numbers (De > 1) the extensional viscosity in a uniaxial extensional ﬂow
scales with the square of the extensibility parameter, L (cf Chilcott and Rallison, 1986

and Bird et al., 1987),

°‘ Lz-l , (2.8)

32

 

 

 

 

2 o M_=0.90-10' T=21°C
s o M. =1.25 -10'
1E+04~1 c1 M.=2-11'10.
: .°° <é> = 16.5 1/s
.4 0. <0) =14.01/S
P _ 9° ..0.
\ 1000E . , . . ,
+ u 2 .° . I 0 I °
P 1 ° °. 0. .
“a . I . . I °
. . . I . (é) 3314.61/8
100's 0 . o I o
i -
.4
1O 1 T f T T f
0.00 050 1 00 1.50 200 250 300
’ t/x,

Figure 2.5 (a) Transient extensional viscosities as measured by ﬁber spinning for the
PIB in PB solutions.

33

 

1E+O4d
1'= 21°C

 

ii - 317,
n - 72
ES
8:
erL

 

 

100 . T . . ....,
1£+05 1E+06

 

T I T r r I

1é+07

Figure 2.5 (b) Apparent Trouton ratios, accounting for the polymer contribution, based
onﬁberspinningmeasurementstalmnaté ~ 15 r‘forthePIBinPBsolutions.

34
where L2 represents the ratio of the mean square end-to—end distance for the fully

extended polymer molecule to its equilibrium value. For a theta system, the mean

squareend-to—end distance of a polymer molecule scales linearly with the molecular
weight.

at M, , (2.9)

The mean square end-to-end distance of the fully extended chain scales with the square

of the molecular weight,

(13> or MS . (2-10)

Thus, by deﬁnition the square of the extensibility of a polymer molecule in a theta

solvent is proportional to its molecular weight,

L2 = (I2)

8 M I (Zell)
2 W

and hence,

.. M, . (2.12)

The ﬁber spinning experiments on the PIB/Decalin solution were performed over
a range of ﬂow rates, ﬁber lengths and draw ratios (see Table 2.2). Figure 2.6 shows
the apparent Trouton ratio plotted against the extension rate for ﬁve different tests.
Here, the average extensional viscosity is divided by the shear viscosity evaluated at a

shearstrainrateofj=ﬁé

2.5 Wm

2.5.1 Permeation apparatus

The experimental apparatus used in this study was similar to that of Chnrielewski
er al. (1990a) with modiﬁcations aimed at minimizing polymer degradation and providing
a means of ﬂow visualization. A set of experiments was initiated by charging the
holding tank with a fresh 8 liter batch of test ﬂuid (see Figure 2.7). Approximately 2
liters of this ﬂuid was then pumped via a peristaltic pump (TAT Engineering, model 110-
43E) into a reservoir. Regulated nitrogen was used to move the ﬂuid from the reservoir,
through the permeation apparatus and back into the holding tank. Pressure transducers
on either side of the cylinder array provided pressure drop measurements, and the ﬂow
rate wasobtained by monitoring the change in liquid level within the reservoir over time.
The same 8 liter charge was used throughout the permeation testing of one cylinder pitch

type. When the cell geometry was changed a fresh batch of ﬂuid was employed.

 

 

 

 

1E+04d
I T = 23 °C
.4
2
it
v
? 10001 '
\ .
It; 1 g
100 I I I r r t r117 1 y ' rrrrr
1 10 100

t’: [3"]

Figure2.6 ApparentTroutonratiosbasedonﬁberspinningdataofthePIBindecalin
solutionwherethevaluesofnaredeterrmnedaty -= 131:.

37

Permeation Apparatus

_L N2

 

Reservoir

Pump

:53 L Cell W

 

 

 

 

 

 

 

 

 

 

Holding Tank

Figure 2.7 Schematic of the permeation apparatus.

38
Two permeation cells are used which differ only in the packing geometries of the

cylinder arrays they contain. Figure 2.8 shows a schematic of the permeation cells. The
test ﬂuid enters the cell through a 2.5 cm I.D. pipe and travels through a 7.6 cm long
transition region before reaching the cylinder array. The cross section of the cell is
rectangular and measures 6.3 cm by 3.8 cm. As shown in the top view on Figure 2.8,
solid wedges ﬁll the comers of the cell to facilitate a smooth transition of the ﬂuid as it
enters and leaves the cell.

The arrays are composed of acrylic circular cylinders, 0.476 cm in diameter and
3.8 cm long. The cylinders are situated so that their axes are perpendicular to the ﬂow,
and they are ﬂush with the cell walls to prevent ﬂuid channeling. Wall and end effects
are known to be negligible (see Chmielewski er al. , 1990a). Cylinder bed lengths in
both cells are approximately 9.3 cm (12 rows) and contain ﬂush mount diaphragm type
pressure transducers (Omega Engineering, model PX102) on either side. A chart
recorder was used to monitor the upstream and downstream pressures over the duration
of each test. Figure 2.9 shows the square and hexagonal pitch geometries. As shown
in Figure 2.9, the cylinders in the square pitch are spaced 0.771 cm from center to
center and in the hexagonal pitch the cylinders form equilateral triangles which are 0.828

cm on a side and 0.717 cm in height.

2.5.2 Flow unsteadiness
Both upstream and downstream pressure traces of the Newtonian ﬂuids showed

no variation with time. This was not the case with the viscoelastic ﬂuids. Small

39

 

 

 

 

 

 

  

      

7/-//////

/

//////-%

I /

              

“a

/ / '/

\

\\\
\\\\\\\\\\\\

   

\\\\\

\\
\\\\\\\\

   

\\\\\\\\\\

\\\\\\\\\\\\\\

\—
\x\\\\\\\\ \\

\\\\
~\\\\\\\\\\\

   

\\\\\\

/////////////////—//////////////
——l 7.6 l'-—9.3“'| 7.6+:T

Side View (cross section)

 

Figure 2.8 Schematic of the permeation cells (all dimensions in cm).

 

 

Figure 2.9 Geometry of the circular cylinders (radius = 0.238 cm) (a) square pitch (b)
hexagonal pitch.

41
amplitude oscillations were recorded by the downstream pressure transducer at a Deborah

number of 1.2 for all of the non-shear thinning elastic ﬂuids. Above this onset, the
pressure oscillations grew in amplitude up to 2 percent of the average pressure with
increasing Deborah number; the Deborah number was controlled by the ﬂow rate. At
no time were pressure ﬂuctuations observed upstream of the cylinder arrays. No
pressure ﬂuctuations were recorded in ﬂow through a blank channel at Deborah numbers
up to 4. This indicates that the unsteadiness resulted from ﬂuid elasticity in ﬂow through
the array of cylinders.

Two examples of pressure traces for the square and hexagonal array of the 2.11
x 10‘ molecular weight solution at De ~ 2.6 are given in Figure 2.10. In both traces
the upstream and downstream measurements are on a different scale. Figure 2.10 (11)
represents the pressure over time of a test run in the hexagonal array. At a ﬂow rate of
81 cm’ls the upstream pressure measured 2.92 x 10’ Pa (42.3 psi), while the downstream
pressure ﬂuctuated around an average value of 0.11 x 10’ Pa (1.6 psi). At this Deborah
number the amplitude of oscillation is approximately 2 percent of the average pressure
and the dominant frequency is approximately 0.2 Hz. A pressure trace at De = 2.6 for
the square array is given in Figure 2.10 (b). Here the ﬂow rate is 113 cm’ls, and the
corresponding upstream pressure is 1.78 x 10’ Pa (27.8 psi). The downstream pressure
oscillates around the mean of 0.17 x 10’ Pa (2.5 psi). Again the amplitude of oscillation
is about 2 percent of the mean pressure, and the dominant frequency is approximately
0.4 Hz. Though this estimate of a. characteristic frequency of the pressure ﬂuctuations
is crude, it still bears pointing out that these frequencies are of the same order of

42

 

 

 

Mair (Pa)

 

 

 

 

 

“in" 2.10 (a) Pressure traces for the 2.11 x 10' molecular weight PIB/PB solution in

hexagonal array at De = 2.6 (Q = 81 cm’ls, upstream P = 29.2 x 10‘ Pa,
dowastr'eaml’ = 1.1 x 10‘ Pa).

43

 

 

 

 

 

 

 

 

Figure 2.10 (b) Pressure tracesforthe2.ll x 10‘ molecular weight PIB/PB solutionin
thesquarearray ”De =2.6(Q =113cm’ls,upstreamP =17.8x10‘Pa,
downstreamP - 1.7x10‘Pa).

44
magnitude as the velocity ﬂuctuations found by McKinley et al. (1991) near the lip of

a 4:1 axisymmetric contraction ﬂow using similar viscoelastic solutions.

2.5.3 Polymer degradation
The permeation apparatus had several design features which had been incorporated

to minimize the severe polymer degradation found by Chmielewski et al. (1990a) - see
Chapter 4. These features included a reduction in the expansion and contraction ratios
at the entrance and exit of the permeation cell, and an increase in the reservoir size,
reducing the frequency in which a batch of ﬂuid is cycled through the cell. Also, lower
molecular weight PIBs were used in this investigation.

In order to assess the extent of polymer degradation the shear viscosity and
storage modulus were measured on ﬂuid samples taken from an 8 liter ﬂuid batch after
it had passed through an array several times. This information was used to determine
the number of runs after which fresh batches of each ﬂuid were required. Very little
degradation was observed for all of the solutions used in this study. This is in contrast

to the extreme amounts of chain scission found when eare was not taken to minimize

degradation (see Chapter 4).

2.5.4 Permeation of Newtonian ﬂuids
The ﬂow resistance of the test ﬂuids passing through the arrays may be

represented by the friction factor, f,

45

a (-AP) Zap e3
f 1 G, (l-e) (2°13)

where p is the ﬂuid density, -AP is the pressure drop, 6 is the void fraction of the bed,
lis the bed length, a is the cylinder radius and G, is the superﬁcial mass ﬂux. The

Reynolds number, Re, is deﬁned by,

ZaGO 1

11 ———(1_£) (2.14)

 

In this study two Newtonian ﬂuids, whose shear viscosities differ by nearly an order of
magnitude, are used to obtain friction factor results for Reynolds numbers ranging from
0.001 to 0.3. Figure 2.11 shows no signiﬁcant difference in the ﬂow resistance of the
Newtonian ﬂuids in the square and hexagonal arrays at the void fraction level of 70
percent. The theoretical calculations of Sangani and Acrivos ( 1982) also predict little
difference in the ﬂow resistance for these two array types at this void fraction level. For
example, the product of the friction factor and Reynolds number, which is inversely
proportional to the transverse permeability, has a theoretically predicted value of 150 for
the square array and 141 for the hexagonal array at a void fraction level of 70 percent.
Figure 2.11 reveals very good agreement between the experimentally measured friction
factors for the two arrays and the theoretically predicted values over the range of

Reynolds numbers studied. This agreement is also evidence of the insigniﬁeance of wall

 

1 E+05

1000

 

A Hexagonal Pitch. Experimental
0 Square Pitch. Experimental
— Theoretical
100 I U U I U I T V U I I I I V I
0.001 0.010 0.100 1 .000

 

1 U UTUIUUI '

Re

Figure 2.11 Friction factor vs. Reynolds number of two Newtonian PB liquids
compared to the theoretical prediction of Sangani and Acrivos (1082) rn square and
hexagonal arrays (closed symbols represent the higher molecular werght PB).

47

and end effects.

2.5.5 Permeation of viscoelastic liquids

The departure of the friction factor from Newtonian behavior for the elastic
liquids is shown in Figures 2.12 (a) and (b). The Re at which the onset of elastic effects
occurs decreases with increasing molecular weight just as it does in the random packed
bed experiments of Kuliclre and Haas (1984). Also, at the same molecular weight the
onset Re is consistently lower for the hexagonal pitch geometry than square pitch.

In order to scale out molecular relaxation time differences related to the PIB
molecular weights, Figures 2.13 (a) and (b) present the relative ﬂuid resistance,

f-Re/(f-Re)”, versus the Deborah number,

V

where v, is the superﬁcial velocity. With this deﬁnition of De, a clear distinction may
be made of the onset of viscoelastic effects between the two array types, 0.80 for the
square pitch array and 0.35 for the hexagonal pitch array. The onset Deborah number,
however, is independent of molecular weight. At sufﬁciently large Deborah numbers the
relative ﬂow resistance becomes independent of the Deborah number. This has also been
observed by Kulicke and Haas (1984) and James and McLaren (1975) for the ﬂow of

polymer solutions through random beds of spheres. The asymptotic value of the

48

 

 

 

 

1E+04,
. , ' . '0...
J -_ . '
J \\\
..eOOOOe
- 10004, "u, ',’.--""'~'.':'°-..,
I o
J
. Square Pitch
10 M,-o.oo-1os
e u, - 1.25 - 10-
100 ° “-72'13'Lve...., . sever-..
0.01 0.10 1.00
Re

Figure 2.12 (a) Friction factor vs. Reynolds number of the PIB in PB solutions ﬂowing
through the square array.

49

1 E+05

 

L lLllll

154-04? \e e .

\ "......

‘”Oe' ' ' '0.
I

I u . .l x
1 000-: 5..
' Hexagonal Pitch

a M, - 0.90 '10'
e M. - 1.25 '10'
100 o 313%.1U1r.11?1 V V r II’WII' ‘ I’ rTlTﬁ

0.00 0.01 0.10 1 .00

L Lilli

IWL Llj

 

 

 

Figure 2.12 (b) Friction factor vs. Reynolds number of the PIB in PB solutions ﬂowing
through the hexagonal array.

 

100.0
‘ o M, - 0.90 - 10' Square Pitch

e M, - 1.25 '10'

u M, - 2.11 no

RAIL

10.01 ...-urn"-
.. ...

a . 00'.

‘35:: "'
.
r

fo Re/(f' Re)"

1.0 . TW

 

 

_l_ Llllll

 

 

 

0.1 . cm”... .
0.1 1.0 10.0

De

ﬁgure 2.13 (a) Flow resistance vs. Deborah number for the PIB in PB solutions
ﬂowingthroughthesquarearray.

51

 

 

 

 

 

100.0 .
3 o u. - 0,90 . 1o- Hexagonal Pitch
1 e M. - 1.25 ' 10'
J D M. II 2.11 " 10'
J ...-I
.. ....” .
I: 10.0? .. ...OO
& E . I. .*.”..00
so.— 1 e) ..
V I
\ .. ”
0 " .ee
0.: 0"
“- 1.01—9-4"——.-”r.
J
0.1 U V U I 1 I T j V
0.1 1.0 10.0
De

ﬁgure 2.13 (b) Flow resistance vs.
ﬂowing through the hexagonal array.

Deborah number for the PIB in PB solutions

52

resistance ratio increases with increasing molecular weight in both arrays, but is
consistently lower in the square array than the hexagonal array at comparable molecular
weights. Figure 2.14 is a plot of the numerieal value of the asymptotic ratio versus
molecular weight of the polymer solute for each array. In both arrays the value of the

ﬂow asymptote seales linearly with the molecular weight,

f°Re
ﬁfe-Elm»: “ M, , (2-15)

just as the apparent Trouton ratio of these solutions does (Figure 2.5 (b)). This is
reasonable beeausc the kinematics of the ﬂow of polymer solutions transverse to cylinder
arrays is dominated by planar extension at high Deborah numbers. In planar extensional
ﬂows, as in uniaxial extensional ﬂows, the extensional viscosity of a FENE model ﬂuid

at high Deborah numbers is proportional to L’, and thus to the molecular weight,

"3.2“;

“La-IO‘M . (2.17)
n-n. '

Hence,

100

53

 

4 0 Square Pitch
A HexagonolPhch

 

 

 

 

0
H
0 a
‘5.
H
§ .§
< .J
m 10-
Q) .
o: D .
‘D .E." i
.2 :1: J
4.1
2
Q)
m J
1 T r r rrrT T V T V ' ' r
1E+05 154-06 1E+O7

Mol. Wt. [g/mol]

I'Igure2.l4 ThehighDeborahnumberﬂowresistanceasymptotesasafunctionof
molecularweightforboththesquareandhexagonalpitcharrays.

54

fits a "ls-2n:
f.ReTN|D.>>1 —'I 5'": ' (2'18)

This relationship between the asymptotic values of the relative ﬂow resistance and the
solution molecular weight comes about as a result of differences in the degree of
extensibility of the polymer molecules. For theta systems the square of the extensibility
is proportional to the molecular weight (Eqs. 2.9-2.11). These results also indicate that
even though ﬁber spinning is an axisymmetric extensional ﬂow, the results are quite
relevant to the planar geometry of ﬂow through cylinder arrays.

Figure 2.15 shows the ﬂow resistance results for the PIB/Decalin solution in both
the square and hexagonal pitch cylinder arrays. Here the Deborah number is deﬁned

With AK‘S'),

I
A ' = —ﬂ‘2—)— . 2.19
1”) wG”(o)) I”? ( )

These data show a plateau at a ﬂow resistance less than 1 for both the square and
hexagonal arrays for Deborah numbers in the range of 0.3 to 1. We were not able to
work at sufﬁciently low Deborah numbers to observe Newtonian behavior. This
reduction in the ﬂow resistance was not observed for the Boger liquids and is a result of
the highly shear thinning nature of the PIB/Decalin solution. As the Deborah number

is increased further, the ﬂow resistance increases, and at De = 3 for the hexagonal array

55

 

 

 

 

 

10.0 .
j 0 Square Pitch
.4 a Hexagonal Pitch
.4
+ a
<5 ‘ f
a: 4‘
C 1.0 He
\ .1 on
a) .
0!. . a
H— '1 AAA A A
4 m D
J
0.1 7 ﬁr r W— 1 r r T r1 f I U I F r I I
0.1 1 .0 10.0
De

Figure 2.15 Flow resistance vs. Deborah number for the PIB in deealin solution ﬂowing
through both the square and hexagonal array.

56

and De = 4.2 for the square array the ﬂow resistance crosses the relative f-Re axis at
1. This enhancement in the ﬂow resistance at high Deborah numbers is a result of
extensional viscosity effects dominating the shear viscosity effects. This is similar to the
ﬁber spinning results where it was found that the averaged extensional viscosity of the

PIB/Deealin solution increased with increasing extension rate.

2.6 Cnnrlusians

The effect of varying polymer extensibility on the dynamics of polymer solutions
ﬂowing transverse to cylinder arrays is studied. The polymer extensibility of each
solution is controlled through the molecular weight of the polymer solute. Average
extensional viscosity measurements afﬁrrn that increasing molecular weight corresponds
to increasing polymer extensibility. The apparent Trouton ratio at a ﬁxed stretch rate
seales linearly with the molecular weight. The pressure drop of non-shear thinning
elastic liquids ﬂowing through both square and hexagonal pitch cylinder arrays at a void
fraction of 70 percent is enhanced at onset Deborah numbers of 0.8 and 0.35
respectively. At high values of Deborah number the relative ﬂow resistances in both
arrays become independent of the Deborah number. This was observed after care was
taken to avoid degradation of the polymer. The magnitude of the high Deborah number
ﬂow resistance asymptote seales linearly with the molecular weight, and hence correlates

well with the apparent Trouton ratios measured.

Chapter 3

THE KINEMATICS OF VISCOUS AND
VISCOELASTIC LIQUID FLOWS WITHIN ARRAYS
OF CIRCULAR CYLINDERS

3.1 Summary

The kinematics and hydrodynamic stability of viscous and viscoelastic liquid ﬂows
transverse to periodic arrays of circular cylinders has been studied at Reynolds numbers
less than 0.5. Both streak photography and laser Doppler velocimetry were used to
observe ﬂow transitions resulting from ﬂuid elasticity in square and hexagonal pitch
arrays at a porosity level of 70 percent. Below an onset Deborah number, the ﬂow of
a non-shear thinning elastic liquid was steady, spatially periodic, and identieal to the
experimentally observed Newtonian kinematics and Stokes ﬂow simulations. LDV
measurements made above the onset Deborah number reveal ﬂow unsteadiness in both
array types. Particle path asymmetry is also observed above the onset Deborah number.
The onset Deborah number corresponds approximately to the onset of elastic effects in

ﬂow resistance measurements found in the previous chapter: 0.70 for the square array

57

58
and 0.25 for the hexagonal array. Results with a shear thinning liquid underline the

elastic origin of the instability observed. Also, these results indicate that any attempt to
predict ﬂow resistance increases must describe the viscoelastic transition to unsteady

ﬂow.

33W

The question addressed in this chapter is whether ﬂuid elasticity, resulting from
the dissolution of small amounts (0.2 wt. 96) of a high molecular weight polymer into
a viscous liquid, affects the ﬂow kinematics within periodic arrays of circular cylinders.
This issue is examined with two experimental techniques - streak photography and laser
Doppler velocimetry (LDV).

The study of how velocity ﬁelds are affected by ﬂuid rheology is important to
both the experimentalist, who uses this information to evaluate ﬂow dynamics, and to the
theorist, who uses the data to evaluate constitutive models and computational methods.
The ﬂow of non-Newtonian liquids past a cylinder provides a good example of this.
Using tracer dyes, Manero and Mena (1981) visualized the slow ﬂow (Re < 0.01) of
shear thinning elastic liquids around circular cylinders. They found that in the range 0.2
< De <1thestreamlinesshiﬁeddownsneamfromﬂresymmeuicpanerncharacteﬁsﬁc
of Newtonian ﬂuids. This was in qualitative agreement with the perturbation calculations
of Mena and Caswell (1974), using an Oldroyd constitutive model. At De ~ 1 there

was no displacement in streamlines, and for De > 1 the streamlines moved upstream of

59
the Newtonian pattern.

In the present study another ﬂow visualization method, streak photography, is
used to qualitatively examine the ﬂow ﬁeld within arrays of cylinders. In this method
small particles in the ﬂuid reﬂect light from a sheet of light illuminating the ﬂow ﬁeld.
Theﬂowisthenphotographedatlongexposuretimessothatthepathsofseveral
particles showupasstreaksontheﬁlm. 'l‘histechniquehasbeenusedextensivelyby
several researchers to visualize ﬂow pattern changes and instabilities, resulting from ﬂuid
elasticity, in axisymmetric entry ﬂows through circular tubes. For example, Nguyen and
Boger (1979) present a series of photographs revealing several ﬂow transitions with
increasing Deborah numbers for non-shear thinning, elastic ﬂuids - Boger ﬂuids - ﬂowing
through a 7.675:l axisymmetric contraction. At sufﬁciently low Deborah numbers (~
0.5),theﬂowpattemwassimilartotheNewtonianﬂuidpattem,havingasmall
secondaryvortexinthecorneroftheupstrearn tube. Astheﬂowrateincreased sodid
thesizeofthevortex, beyondwhatisseenwithNewtonianliquids. AtDe ~ 3the
vortex continued to grow and became asymmetric. At De ~ 6 the asymmetric vortex
began to rotate around the tube. Finally, at De ~ 15 the ﬂow beeame chaotic.

In addition to ﬂow visualimtion, laser Doppler velocimetry (LDV) has been used
to obtain quantitative local velocity measurements within cylinder arrays. LDV is an
experimental technique which measures point velocities within a ﬂuid ﬂow by detecting
the Doppler shift of light scattered from solid particles in the ﬂuid. Presumably these
particles are small enough so that their velocity corresponds to the loeal ﬂuid velocity.

Since this technique’s introduction by Yeh and Cummins (1964), LDV has been

60
widely used on both laminar and turbulent ﬂows of gases and Newtonian liquids. The

use of LDV on ﬂows of viscoelastic liquids, however, has not been as extensive. An
LDV study particularly relevant to this investigation has been made by Lawler at al.
(1986). They examined the velocity ﬁeld and ﬂow transitions of an elastic, non-shear
thinning polyisobutylene solution (a Boger ﬂuid) in a 4:1 axisymmetric contraction ﬂow.
For De < 0.8 the ﬂow was steady and identical to the predicted Newtonian ﬂow ﬁeld,
but at De = 0.8 the ﬂow became time periodic with a ﬂuctuating tangential velocity
component. At De = 1.2 the ﬂow again beeame time independent, but was no longer
identical to the Newtonian velocity ﬁeld. It is interesting that the ﬂow transition at De
= 0.8 is lower than the onset Deborah number for the appearance of large corner
vortices, indieating that the LDV technique is sensitive to ﬂow transitions, particularly
temporal transitions, which otherwise could not be detected by more conventional ﬂow
visualization techniques.

The work presented here examines the kinematics of viscous and viscoelastic
liquids in periodic arrays of circular cylinders. The ﬂow ﬁeld is explored by both streak
photography and by laser Doppler velocimetry. LDV provides point velocity
measurements of both viscous and viscoelastic liquids in cylinder arrays for comparison
with numerieal simulations, while streak photography provides global comparisons of
ﬂows at Deborah number of 0(1) with the Newtonian ﬂow ﬁeld. In both the square and
hexagonal arrays a ﬂow transition is observed for the viscoelastic liquid near the Deborah

number where onset of elastic effects is observed in the ﬂow resistance data.

61
3.3 W

3.3.1 Basic principles

A dual beam, single component LDV system is used in this study, and is shown
schematieally in Figure 3.1. Plane polarized light, emitted from a 35 mW He-Ne laser
at a wavelength of 632.8 nm, is split into two beams of equal intensity and focused by
a transmitting lens. The receiving optics and photomultiplier are set up in the forward
seatter mode, directly in line with the transmitting lens. It is in this direction that the
intensity of seattercd light, resulting from small particles moving through the beam
intersection, is the greatest. Velocity measurements are made in the ellipsoid, known as
the “measuring volume", formed by the beam intersection and it is in this location,
between the transmitting and receiving lenses, that the permeability cells were placed.

The simplest and most widely used explanation of the operation of the dual beam
LDV system is based on the fringe model. This model avoids reference to the Doppler
shift effect and yet provides many correct results in terms of the velocity measurement.
It breaks down, however, in the calculation of predicted signal intensity of light seattered
from the measuring volume.

The fringe model is based on the interference of intersecting waves at the

measuring volume, producing a fringe pattern with spacing, df,

62

 

.:.m use... 5
coca—.02 0.8 3:28an8 82%.." 3.25029, 33qu .53— .8 3882.8 ﬁn Bear—

 

 

 

 

 

 

 

 

 

—F-

 

 

 

 

63

Table 3.1 Components of the laser Doppler velocimeter system shown in Figure 3.1.

 

 

 

 

 

 

 

 

 

 

 

Component TSI Model Description
1 9126-255 Laser power supply
2 9126-105A 35 mw lie-Ne laser, A 632.8 nm
3 9115-2 Beam Splitter, 50 mm separation
4 9118 Transmitting lens, f = 250 mm
5 9165 Photomultiplier power supply
6 9118 Receiving lens, f= 250mm
‘7 9160A Photomultiplier
8 9126, 9121 Optical rails
9 1980 Signal Processor, 100 MHz clock
10 465 M Portable Oscilloscope
11 ~—-- Apple IIe computer

 

 

64

A

Zainht) (3.1)

f.

where x is the wavelength of the laser light and x is the half angle of the intersecting
beams (x = 5.71 °). In the present system the fringe spacing is 3.18 um, resulting in 64
fringes within the measuring volume. These fringes are set in planes perpendicular to
the plane in which the beams lie, and run parallel to the line bisecting the angle formed
by the beam crossing. Figure 3.2 shows a schematic of the ellipsoid and fringe pattern
at the beam intersection. Here the beams cross in the x-z plane and the z axis lies along
the bisector of the angle formed by the beams.

The velocity measurement of ﬂuid ﬂowing through the measuring volume relies
on small solid particles, traveling with the ﬂow, to scatter light. As a particle moves
through the fringe pattern, past the light and dark hands, it reﬂects light with an
oscillating intensity. This reﬂected light is picked up by the photomultiplier and
converted into an electronic signal. The signal, shown in Figure 3.3, contains a mean
low frequency component known as the pedestal and a sinusoidal component which
oscillates at the Doppler frequency (ran). During processing of the signal, the pedestal

is usually removed, leaving an which is proportional to the velocity,

v .Jﬂ . (3.2)

D d!

 

 

 

 

 

Figure 3.2 Measuring volume and fringe pattern formed by the beam intersection. The
beams intersect in the x~z plane and the z axis follows the bisector of the angle of beam
crossing (Dabir, 1983).

Intensity —' ‘— “ ' ‘f"‘

 

 

3:331; a.“ ll ‘ ‘ l

”H - new:
/ ”HAHN" “l" H. / ,.

 

 

Time

Figure 3.3 A typical signal from the photomultiplier (Dabir, 1983).

67
Hence, calculated ﬂuid velocity is based solely on the measured Doppler frequency, the

wavelength of laser light and the half angle of the beam crossing. The dual beam system
measures only velocity components perpendicular to the plane of the fringes (the x
direction in Figure 3.2) since the other components represent particle motion within a
fringe plane, producing no oscillating signal.

The absolute value in Eq. 3.2 indieates this technique’s inability to determine
ﬂow direction. A particle traveling through either end of the measuring volume at the
same speed will result in an identical signal. In ﬂows where the loeal direction of ﬂuid
motion is unknown a technique eallcd frequency shifting is employed. Instead of the
measuring volume containing a stationary interference pattern, the frequency of one of
the transmitted beams is shifted by 0,, causing a moving wave-like fringe pattern. Thus,
particles in the measuring volume moving in the same direction as the fringes will result
in experimentally measured frequencies less than u,, while particles moving in the
opposite direction will result in frequencies greater than 0,.

3.3.2 Experimental LDV system

Figure 3.1 shows the dual beam Thermo-Systems Inc. (TSI) LDV system used
in this study. Each of the individual components is listed and described in Table 3.1.
All experiments are performed in the forward scatter mode and frequency shifting is not
employed. The light source is a 35 mw He-Ne laser, emitting light at a wavelength of
632.8 nm. The major optical components consist of a beam splitter, which divides the

incident beam into two equal intensity beams separated by 50 mm, and transmitting and

68
receiving lenses, each having a focal length of 250 mm. This lens system produces a

beam halfangle of 5.l7° and an ellipsoid measuring volume with dimensions: A2 = 1.9
mm, Ax = 0.18 mm and Ay = 0.18 mm (see Figure 3.2). As shown earlier by Eq.
3.1, the fringe spacing is 3.18 um, resulting in 64 fringes across the measuring volume.
Also, the industrial grade solvents used in this study are sufﬁciently contaminated with
dust and other forms of dirt that seeding the flow is unnecessary. The velocities
measured in this study typically are in the range of 0.5 to 5 cm/s, which correspond to
Doppler frequencies of approximately 1.6 to 16 kHz.

The processing of the Doppler signal begins at the photomultiplier which picks
up photons from the receiving optics and converts them into a voltage signal. The
voltage signal is then sent to a TSI model 1980 signal processor where it passes through
an input conditioner and a timer. The input conditioner ampliﬁes and ﬁlters the signal
with 1 kHz and 10 kHz (or 100 kHz) high and low pass ﬁlters. The conditioner also
contains a Schmitt trigger. If the signal amplitude is greater than 50 mV the trigger is
activated, converting the sinusoidal wave into a square wave. Otherwise, the output is
not updated. The timer’s function is to measure the length of the envelope containing
N cycles from the Schmitt trigger; the numbers of cycles per burst, N, is set externally
andforthisstudyisS. Thetimeralsomeasuresthelength ofanenvelochlchcles
long. If the average signal frequency of the ﬁrst N/2 cycles is not within 5 percent of
that of the N cycles then the data point is rejected and the system is reset without

updating the output. Otherwise, the frequency is latched to output.

3.3.3 LDV meamrernent difficulties

The LDV measurements, made on both the viscous and viscoelastic liquids within
the cylinder arrays, suffered from low data rates. Data rates of approximately 1 to 5 Hz
wereobtained. ThisisincomparisontoratesoflOOtoSOOszhicharenecessaryin
order to obtain temporal information about the velocity. As a result, sampling was
limited to the ”handshake" mode of data collection. Typically, at data rates of 100 Hz,
data can be collected every 100 ms with conﬁdence that every data point is independent.
\Vrth a data rate of 1 Hz and sampling every 100 ms, however, every tenth data point
is independent. The other nine represent the same point because the previous signal will
remain in memory until a new point replaces it. In the handshake mode of operation data
points are sampled only as frequently as they arrive. Thus, making accurate velocity
measurements of steady ﬂows was not a problem. Obtaining temporal information from
unsteady ﬂows, however, was not possible as a result of the low data rates.

During this investigation, sets of 256 data points were sampled at each location
intheﬂowﬁeld;atdataratesoflHzthisprocesstookmorethan4minutesperset.
A TSI data reduction program, running on an Apple IIe, collected the frequencies from
the signal processor, calculated the corresponding velocities and presented a statistical
evaluation of the data which included the probability distribution function, the mean
velocity and standard deviation. Collecting data this way worked well since the standard
deviations of the data sets were less than 1 percent of the mean values for the stable
ﬂows.

Theprecisecauseofthelowdataratesisunknown. Whatisknownisthatthe

70
low rates were not a result of infrequent Doppler bursts, but rather small signal

amplitudes. The amplitude of the electronic signal received from the photomultiplier was
usually lower than the 50 mV necessary to activate the Schmitt trigger. Thus, the timer
infrequently received data to compare and latch to output. This occurred even though
the signal gain was set at maximum. The cause of this low output amplitude was initially
thought to be the result of a faulty photomultiplier. However, a TSI inspection has
revealed that this is not the case. Furthermore, the entire LDV system was tested by
measuringthevelocityofwaterstirredinaglass beaker. GoodDopplersignals were
obtained with data rates of 100 Hz. Seeding the ﬂow with spherical 5 micron Nylon
particleswasalsotested. ThisappearcdtoslightlyincreasethefrequencyofDoppler
bursts, but had no effect on the signal amplitude.

The most likely cause of the low signal intensities is optical inhomogeneities in
the Plexiglas windows of the permeability cells. Because the cylinders are tightly press
ﬁt into the windows, much residual stress remains in the material. This is conﬁrmed by
the asymmetric stress patterns observed in the Plexiglas when they are examined between
two polarizing lenses. These patterns are the result of local index of refraction
differences caused by residual stresses. Thus, light exiting the permeability cell is
scattered by the Plexiglas window, reducing the intensity of light picked up by the

receiving optics.

71
3.4 Maternal

3.4.1 Test ﬂuids

Three test ﬂuids were used in this study: a Newtonian ﬂuid, a non-shear thinning,
elastic ﬂuid (a Boger ﬂuid) and an elastic shear thinning ﬂuid. The Newtonian liquid is
a pure polybutene (PB), Amoco grade H25, whose material and rheological properties
can be found in the preceding chapter. The Boger liquid is a polyisobutylene! kerosene!
polybutene mixture, having the following composition by weight: 0.25 % PIB/ 7 96
kerosene! 92.75 % PB H25. The PIB is a 4 to 6 million molecular weight polymer
purchased from Aldrich Chemical Co. This solution is a Boger ﬂuid similar to the M1
standard (Sridhar, 1990) with the shear viscosity nearly constant over shear rates up to
10 s" (see Figure 3.4). The shear viscosity, shown in Figure 3.4, is 2.8 Pa-s at 25 °C
and the relaxation time, x“ calculated from the quadratic region of the storage modulus
vs. frequency curve (see Figure 3.5) is 0.86 s. The shear thinning liquid is a 2 wt.%
PIB in decalin solution whose rheological properties are presented in the preceding
chapter.

3.4.2 Apparatus
Theﬂowloopandperrneabilitytestcellshavebeendescribedindetailinthe

previous chapter, and only the important features of the apparatus are discussed here.
AschematicoftheapparatusisshowninFigureZJ. Testﬂuidispumpedviaa
peristaltic pump from the holding tank into the liquid reservoir. From here the test ﬂuid

 

 

 

 

10
i T = 25 “c

.—

‘1’

O

2:.

c e e . . . . 0
1 r I t at l "I r v ...
0.01 0 10 1.00 10.00

1 [s41

Figure 3.4 Steady shear viscosity of a 0.25 % PIB/PB solution at 25 °C.

73

 

 

 

1000.00-
0 0' T .. 25 '0
. GUI
100.00 . . ' "
r; e ° . . . e 0
O. 10.00 . g ' . . e .
L—l . O . . 0
ED . O . . . e . .
_- 1.00 . . - . . -
0 . e ' . e '
0.10 . o .
0001 I V V 1 V 1 VIII I U err
0.1 10 10.0 100.0

Figure 3.5 Linear viscoelastic properties of a 0.25 96 PIB/PB solution at 25 °C.

 

74
either ﬂows freely or is forced by regulated nitrogen pressure through the cells. Flow

rate measurements are made by weighing ﬂuid samples collected over time or by
monitoring the reservoir liquid level with time. Pressure transducers on either side of
thecylinderarray measurethemeanpressuredropacross thebedasshowninFigure
2.7.

The kinematics of ﬂow transverse to a square array and to a hexagonal array of
circularcylindersis studied. Each arrayhasavoid fractionof70percentandis
composed of acrylic cylinders (radius = 0.238 cm) arranged as shown in Figure 2.9.
Localized velocity measurements within the cylinder arrays are made by LDV and have
been discussed above.

3.4.3 Streakphotography

Theﬂowpattemsaremappedbystreakphotography. Thisisdonebypassing
a light beam (from the 35 mW laser used in the LDV experiments) through a cylindrical
lenstoformathinsheetoflight. Thissheetisthenreﬂectedoffasurfacecoatcdmirror
andpassedthroughawindowinthetopofthepermeabilitycell,illuminatingtheﬂow
ﬁeld. Thelightsheet liesperpendiculartotheaxes ofthecylinders andintheplaneof
the2dimensional ﬂow. Picturestakenwitha35 mm camera, usingatime exposure of
1.5 seconds,capturethepathof 50 pm silicon carbide particlesseeded(0.033 gramsper

liter ﬂuid) in the ﬂow.

75
3.5 We

3.5.1 Stokes flow simulation

Finite element simulations have been carried out for two dimensional Stokes ﬂow
inthedomainsshowninFigure3.6. Thesedomainsrepresentrepeatunitsofthe
periodic square and hexagonal cylinder arrays used in the experimental section of this
chapter. Both the x and y coordinates are sealed with the cylinder radius and the velocity
with the superﬁcial velocity. Over the entire boundary of both domains the y component
of the velocity is zero. This is the result of no-slip at cylinder surfaces and symmetry
along all other boundaries. The x component of the velocity is set to zero only at the
cylinder surfaces. Otherwise, it is speciﬁed at the upstream portion of the domain
boundary (x = -1.62 and x = -3.01 for the square and hexagonal pitch geometries
respectively). The periodicity of both geometries requires that the x velocity component
along the downstream domain boundary be equivalent to the upstream boundary. This
condition is satisﬁed by iterating; the boundary velocities upstream are replaced by
calculated downstream velocities until convergence is attained. The results of the
simulations are shown in Figure 3.7. As expected for linear ﬂuids with periodic
boundary conditions, the streamlines in both the square and hexagonal arrays are
symmetric around the cylinders.

3.5.2 Newtonian liquid flow vbualization

Theﬂowpanernsforthepmepolybumneinboththesquareandhexagonalanays

76

 

 

 

 

 

 

 

 

 

—— 1.62
\ ’-
l l l I I "
162 1 0 1 182
(a)
y l
1.74
_' 0.74
g ’-
l l l I l x
-3.01 -1 0 1 3.01
(b)

Figure 3.6 Geometric domains used for the Stokes ﬂow simulations (a) square array and
(b) hexagonal array.

 

(a)

 

(b)

1“igure 3.7 Streamline output from the Stokes ﬂow simulations (a) square array and (b)
hexagonal array.

78
at Re = 0.027 and Re = 0.013 respectively, are shown in Figure 3.8. The ﬂuid is

travelling from right to left. The particle paths in the photographs are symmetric,
reﬂecting the periodicity of the arrays, and are identical to the streamline calculations in
Figure 3.7. In the square array no streaklines are visible in the space between cylinder
stagnation points because of extremely low velocities (compared to the bulk ﬂow)
between rows of cylinders. Streak photographs for Reynolds numbers up to 0.5 in both
arrays have also been taken but are not presented here. These photos show the same
patterns as those in Figure 3.8. Also, ﬂow resistance measurements for these tests are
in good agreement with Darcy’s law (see Chapter 2), and at no time did the ﬂow exhibit

any instability.

3-5.3 Newtonian liquid LDV measurerrrents
The results of LDV measurements taken along lines of geometric symmetry, y =
l .62 and y = 1.74 for the square and hexagonal arrays respectively, are shownin Figure
3-9. Along these lines the geometric symmetry requires that the y component of the
Velooity vanish at steady state. Again, good agreement is found between the numerical

Simulations and the velocity measurements.

3-6 WM

345.1 Flow visualization on the PIB/PB liquid

Figures 3.10 (a)-(e) and 3.11 (aHd) show a series of streak photographs, each

... Hum ‘t ‘,.‘

 

 

 

L

l“igure 3.8 Streak photographs for a Newtonian polybutene liquid in the cylinder arrays
at a void fraction of 70 percent (ﬂow from right to left) (a) square array, Re = 0.027
and (b) hexagonal array, Re = 0.013.

 

(b)

4.5

 

-— Numerical
0 Re - 0.001
4.05 0 Re - 0.011

 

 

 

 

3.5~
< 30~
>
J
2.5-
2.0-
J
1.5 . . . ,

V I 1 r V l U T I
-1.75 -1.25 -.75 -.25 0.25 0.75 1.25 1.75

Figure 3.9 (a) LDV measurements and the Stokes ﬂow prediction for the inelastic
Newtonian ﬂuid in the square array along y = 1.62.

I

81

 

 

 

 

5.
~ — Numerical
" 0 Re - 0.011
0 Re - 0.025
4.-
" 3” O
. . ' l
.3.“ e
>° '
\ ..
> d .
2.—
'1 e
1.~
O f r l v v I e r y 1 '
-2 0 -1.0 0 0 1 0 2.0

Figure 3.9 (b) LDV measurements and the Stokes ﬂow prediction for the inelastic
Newtonian ﬂuid in the hexagonal array along y = 1.74.

82

at a consecutively higher De, of the PBIPIB solution ﬂowing through the square and
hexagonal arrays. In all photographs the ﬂuid is travelling from right to left and the
Reynolds number is less than 0.5. The relative ﬂow resistances, corresponding to the
ﬂow situation in each photograph, are given in Figure 3.12. The ﬂow resistance data
are consistent with those presented in the previous chapter, showing an elastic onset at
De = 0.70 in the square array and at De = 0.25 in the hexagonal array. Sufﬁciently
large Deborah numbers were not attained in either array to observe the asymptotic high
Deborah number limit of ﬂow resistance for this high molecular weight PIB.

The photos in Figures 3.10 (a) and 3.11 (a) were taken of ﬂows at De = 0.16
and De = 0.06 respectively. As seen in Figure 3.12 these Deborah numbers are below
the onset values. Accordingly, the particle paths in these photos are symmetric and
match both those found with the Newtonian ﬂuids and those calculated in the computer
simulation. As the Deborah number increases past the onset values, the ﬂow in both
arrays go through a transition from a steady Newtonian ﬂow to an unsteady ﬂow. This
has been observed in both arrays through the downstream pressure ﬂuctuations discussed
in the previous chapter and by LDV measurements discussed in the next section. The
streak photographs presented here also capture this unsteadiness.

In the case of the square array, both ﬂow unsteadiness and asymmetry are
observed. The progression of the asymmetry with Deborah number can be observed in
the photographs shown in Figures 3.10 (b)-(1). This asymmetry is characterized by
particles which follow ﬂow paths winding randomly between cylinders and crossing

through lines of geometric symmetry. This becomes more apparent as the Deborah

83

(a)

24‘

l‘.

. i"

. 12*“ (i

I
3,;
‘4

l
e

 

\,.J

 

'J

1’.

‘ 11‘

..4
1"
is!
I“
’S

...,Q
at

8
a

 

(b)

Figure 3.10 Streak photographs for the viscoelastic liquid in the square array having a
void fraction of 70 percent (ﬂow from right to left) (a) Re = 0.039, De = 0.16 (b) Re
= 0.19, De = 0.80.

 

84

 

 

00

Figure 3.10 (c) Re = 0.25, De = 1.09 (d) Re = 0.36, De = 1.44.

85

 

 

(0

Figure 3.10 (e) Re = 0.48, De = 1.91 (f) Re = 0.48, De = 1.91.

86

 

 

 

 

(b)

Figure 3.11 Streak photographs for the viscoelastic liquid in the hexagonal array having
a void fraction of 70 percent (ﬂow from right to left) (a) Re = 0.015, De = 0.059 (b)
Re = 0.069, De = 0.28.

 

87

5..
.../4r. ........./..../.,.

:3

...». 2., C .53.

e?

s

134
... A

~
\\\\.n./.1..;_1.. mgr

0
....\.
.. ...; . an a
__ .

 

(e)

is. ...... ... ...»... ... ,

e; . _
.3 . near/11%. \

.Em .
i... .. e

..3 .

 

 

(d)

= 0.14, Dc = 0.56.

0.43 (d) Re

,De=

0.11

Figure 3.11 (c) Re

 

 

 

 

 

100.0

. 0 Square Pitch

3 a Hexagonal Pitch

.1 . A

A
1 . - ‘

’5 ° °a .
o: : t °
...... ‘1 a
V .1
\
Q.) “ a
0:
4'— 1.0: ‘ as

1

J

.J

0.1 U T I I'UIVI I I TrUIIUI I I ﬁU'III
0 01 0.10 1.00 10.00
De

Figure 3.12 Relative flow resistance of a 0.25 % PIB/PB elastic liquid in both square
and hexagonal cylinder arrays.

89
number is increased. The ﬂow unsteadiness is apparent by comparing streak lines in

photos taken at the same Deborah during two different experiments. From one snap shot
to the next the particle paths differ - even at the same Deborah number! For example,
Figures 3.10 (e) and (f) are photographs taken oftwo separate test runs, both at De =
1.91 and both at the same location in the ﬂow ﬁeld. Besides noting the asymmetry of
the streak lines, one can also observe identical regions in the ﬂow ﬁeld where the particle
paths are completely different from each other, indicating ﬂow unsteadiness.

Particle path asymmetries and ﬂow unsteadiness are not as easily observed in the
photographs taken here of the high Deborah number ﬂows with the hexagonal array (see
Figure 3.11 (b)-(d). Close examination of the photos reveal many areas where the
particle paths cross each other. This is most apparent in the highest Deborah number
ﬂow shown in Figure 3.11 (d). Particle path crossing is evidence of ﬂow unsteadiness

in the hexagonal array.

3.6.2 LDV measurements on the PIB/PB liquid

Laser Doppler velocimetry measurements of the PIB/PB liquids made along y =
1.62 and y = 1.74, in the square and hexagonal arrays respectively, conﬁrm that for
Deborah numbers below the onset values the kinematics of the ﬂow are identical to those
of Stolm ﬂow (see Figure 3.13). This is no longer true when the Deborah numbers
exceed the onset values.

The LDV measurements made in both arrays at Deborah numbers above the onset

values resulted in Doppler signals which appeared extremely ”noisy" . The Doppler burst

90

4.5

 

. — Numerical (Newtonian)
0 De - 0.63

4.0—

3.5-

q

3.0~

v/v.

2.5-

d

2.0-

 

 

 

 

1.5 1 l U l I I I I

. . . I .
-1.75 -1.25 -.75 -.25 0.25 0.75 1.25 1.75

Figure 3.13 (a) LDV measurements and the Stokes flow prediction for the 0.25 %
PIB/PB elastic liquid in the square array along y = 1.62.

91

 

 

 

 

5.
4 — Numerical (Newtonian)
. 0 De - 0.20
4.4
3.4 e
>0 d ' e
\ d
> .
2.- . ' ' ‘.
d O
1.4
J
o U ﬂ r U U U 1 U U U I U U
—2 0 0. 1 2

Figure 3.13 (b) LDV measurements and the Stokes ﬂow prediction for the 0.251 %
PIB/PB elastic liquid in the hexagonal array along y = 1.74.

92
wasnotcomposedofasinglefrequencyasillustratedinFigure3.3butcontained
multiple frequencies. A local velocity unsteadiness will produce this result. Multiple
scattering particles within the measuring volume, each moving at different velocities, will
result in a Doppler signal containing several frequencies, as illustrated in Figure 3.14
(b). Figures 3.14 (a) and (b) show the probability distribution of 1024 velocity data
pointstakenatx = -1.58andy =1.74inthehexagonalarrayattwodifferentDeborah
numbers, De = 0.085 and De = 0.36 respectively. The extremely narrow velocity
distribution shown in Figure 3.14 (a) indicates the steadiness of the ﬂow below the onset
Deborah number. Above the onset Deborah number, the velocity distribution is very
broad (Figure 3.14 (h)), demonstrating the ﬂow unsteadiness. Similar trends in the
velocity probability distributions below and above the onset Deborah number were seen
in the square array also. As a result of the handshake mode of data collecting it is
impossible to determine whether the ﬂuctuating velocity measurements in either array

were time periodic.

3.6.3 LDV measurements on the PIB/decalin liquid

Figure 3.15 (a) shows velocity data in the square array along the symmetry line
y = 1.62 for De = 0.69 and De = 0.82. At these Deborah numbers the ﬂow is steady,
and the shear thinning properties of the ﬂuid dominate the ﬂow dynamics, resulting in
relative ﬂow resistances below the Newtonian value (see Figure 2.15). The velocity
proﬁles of Figure 3.15 (a) show higher velocities than the predicted Newtonian values

along y =1.62. This is consistent with thereduced pressure dropobserved atthese

93

 

 

 

 

 

 

 

 

 

 

 

 

 

 

   

 

 

s. “WEE-83f UELQCITY PDF ”.3
4
3
2
1
6 LL a

9 .64 rt/S .39

(a)
... FS=EE-t33 l.,.lE|_t:Ir:;IT*¢‘ F’DF ”3,5:
4
3
2 l in
l - ‘ r
illlll‘i’iil 1"" l N

a r

B

 

.84 1‘1 /8 .08

(b)

Figure 3.14 Probability distribution function of the elastic ﬂuid’s velocities measured
in the hexagonal array at x = -l.58 and y = 1.74 (a) De = 0.085 (b) De = 0.36.

94

5.5

 

.. -— Numerical (Newtonian)
5.0 .1 0 De - 0.69
4 0 De I- 0.82

v/v.
U
01
l

 

2.0-4

 

 

105 U I U l' I T" I
-1.75 -1.25 -.75 -.25 0.25 0.75 1.25 1.75

 

Figure 3.15 (a) LDV measurements and the Stokes ﬂow prediction for the PIB/Decalin
liquid in the square array along y = 1.62.

95

 

 

 

5.0
. — Numerical (Newtonian)
0 De - 0.63
" 0 De I 0.94
4.0-
.. . . .
3.0~ . - ' -
>0 " I
I
> . -
2.0- . '
1.0-
0.0 U U U ' U U U I U U U l U U U .
~2.0 -1.0 0.0 1.0 2.0

 

Figure 3.15 (b) LDV measurements and the Stokes ﬂow prediction for the PIB/Decalin
liquid in the hexagonal array along y = 1.74.

96
Deborah numbers. As the Deborah number is increased above 1.0 the extensional

properties of the ﬂuid begin to dominate and the relative ﬂow resistance increases. As
this occurs, a velocity transition similar to that found with the Boger liquid, is observed.
For De > 1, the LDV signal becomes “noisy”, indicating a transition from a steady to
an unsteady ﬂow.

Figure 3.15 (1)) presents the velocity proﬁles along the line connecting two
stagnation points (y = 1.74) in the hexagonal array for De = 0.63 and De = 0.94.
Along this symmetry line the velocities agree with the Stokes ﬂow prediction even though
at these Deborah numbers the ﬂow resistance curve of Figure 2.15 reveals a large shear
thinning effect. As a result of the stagnation points on this particular symmetry line, the
ﬂow along this line is dominated by elongational properties. Shear thinning effects must
occur in the narrow gaps that are not aligned in the x-direction. As the Deborah number
increases above 1, extensional effects begin to dominate the ﬂow and the relative ﬂow
resistance begins to increase. At this point the ﬂow becomes unsteady.

The signiﬁcance of these results is that at ﬂow conditions where shear thinning
effects dominate the ﬂow is steady, even at Deborah numbers as high as 1. It is not until
the extensional effects begin to dominate the ﬂow at higher Deborah numbers that the

ﬂow becomes unsteady.

3.6.4 Discussion of flow unsteadiness
That the onset of ﬂow unsteadiness occurs at the same Deborah number as the

onset of the excess pressure drop suggests that the nature of the two phenomena are the

97
same. In the previous chapter it has been shown that increases in the relative ﬂow

resistance are consistent with and can be correlated by extensional viscosity increases
occurring in localized regions within the cylinder arrays. This is evident in the
hexagonal array. Stokes ﬂow calculations show that along streamlines connecting
stagnation points extension rates of up to 3.5 times the nominal shear rate are obtained.
The large extension rates coupled with the high ﬂuid residence time near the stagnation
points cause macromolecules in these areas to elongate several fold and even break. This
is supported not only by the chain extension calculations of Chmielewski et al. (1990),
but also by the birefringence data of Cressely and Hocquart (1980). Cressely and
Hocquart (1980) studied the ﬂow of polymer solutions around a circular cylinder and
found very localized birefringence along the streamline emanating from the stagnation
point on the downstream side of the cylinder. In the case of the square array the region
of the highest extension rates is localized along the symmetry line running between
cylinder rows (y = 1.62).

At comparable Deborah numbers, the maximum extension rate in the square array
is lower than that in the hexagonal array, resulting in differences in the onset Deborah
number for the two arrays, 0.25 and 0.70 respectively. This difference in onset Deborah
number can be accounted for quantitatively by redeﬁning the Deborah number in terms

of the maximum array extension rate instead of the nominal strain rate, 4,,

De‘ = 1.1-(«7.) , (3.3)

98
where a is 1.5 for the square array and 3.5 for the hexagonal array based on Stokes ﬂow

calculations. Thus, the onset of elastic effects occurs at De‘ ~ 1 for both arrays. This
value corresponds to the Deborah number where macromolecules in extensional ﬂows
undergo the transition from a coiled to an elongated state - the coil to stretch transition
(De ~ 1). Hence, at De' ~ 1 macromolecules in the polymer solutions become
elongated in localized regions within the cylinder arrays, resulting in an excess pressure
drop across the arrays.

The complex ﬂuid dynamics resulting from high extensional stresses generated
in the extensional ﬂow regions may result in the observed ﬂow unsteadiness. An
indication of the apparent Trouton ratio attained in the cylinder arrays can be obtained
from the ﬁber spinning experiments of the previous chapter. The ratio of the average
extensional viscosity to the shear viscosity of a 2.11 million molecular weight pm in PB
solution is 1500 at De = 3.6. Also, data on a PIB in decalin solution show that this
ratio increases with increasing Deborah number.

The onset of ﬂow unsteadiness, occurring at the same Deborah number where
localized regions of large stress appear is consistent with the observations reported by
Ambari et al. (1984) on laminar ﬂow around a single cylinder. They used an
electrochemical technique to study the mass transfer from a circular cylinder in dilute
polyethylene oxide solutions. At an onset Deborah number of approximately 3 they
observedalargedecreaseinthemasstransferratewithrespecttotheNewtonianvalue.
This decrease was accompanied by an onset of ﬂuctuations of the limiting diffusion

current. TheRMSvaluesoftheseﬂuctuationsincreasedwithincreasingchorah

99
numbers, but eventually reached a plateau. Ambari er al. ( 1984) attributed these

ﬂuctuations to the high extensional viscosities localized near the upstream cylinder
stagnation point.

Similar to the observations of Ambari et al. (1984), LDV measurements taken at
a single point in the ﬂow ﬁeld demonstrate an increase in the ”degree" of ﬂow
unsteadiness with increasing Deborah number. The degree of unsteadiness is quantiﬁed
by the ratio of the standard deviation of a set of velocities measured at a point to the
mean value. As the Deborah number of the ﬂow is increased the standard deviation ratio
increases. This is shown in Figure 3.16 which represents LDV data taken in the
hexagonal array at a point approximately one half cylinder radius behind a stagnation
point (x = -1.58 and y = 1.74). Figure 3.16 shows that at De ~ 0.2 we ~ 1) there
is a large jump in the degree of unsteadiness, conﬁrming that this Deborah number
indeed represents a transition point for viscoelastic ﬂow in the hexagonal array. Unlike
the observations of Ambari er al. (1984), the standard deviation of the ﬂuctuations does

not level off for the Deborah number range shown in Figure 3.16 .

3.7 Concussion:

The kinematics of viscous and viscoelastic liquids ﬂowing through square and
hexagonal cylinder arrays has been studied. Both streak photographs and LDV
measurements indicate that below an onset Deborah number, De' ~ 1, the ﬂow of non-

shear thinning elastic ﬂuids is identical to the Stoles ﬂow ﬁeld. However, above the

100

50

 

%

40-1 .

Std Dev/ Mean Vel
N
O
1
.

 

 

 

0.01 0.10 1.00

Figure 3.16 Growth of the ﬂow instability in the hexagonal array for the 0.25 %
PIB/PB elastic liquid as measured by LDV at x = -l.58 and y = 1.74.

101
onset Deborah number the ﬂow becomes unstable. The results with the shear thinning

liquids underline the elastic origin of the instability observed.

The steady Newtonian ﬂow changes to an unsteady and spatially aperiodic ﬂow
ﬁeld at an onset Deborah number where elongational ﬂow dominates. This onset
Deborahnumbercorresponds totheonsetofﬂuidelasticity effectson thepressure drop
across the array. Thus, any attempt to predict ﬂow resistance increases must describe

the viscoelastic transition to unsteady ﬂow.

Chapter 4

THE DEGRADATION OF POLYMER SOLUTIONS
FLOWING THROUGH ARRAYS OF CIRCULAR
CYLINDERS

This chapter was published in the Journal of Non-Newtonian Fluid Mechanics 35, 309-

325 (1990), with eo-authors C.A. Petty and K. Jayaraman.

4.1 Sumarx

The ﬂow of a dilute solution of polyisobutylene in polybutene transverse to
unidirectional arrays of cylinders has been investigated at Reynolds numbers less than
0.1. No different arrays were used - a triangular pitch array and a rectangular pitch
array. Both arrays have a porosity of 0.704, the same bed length and comprise identical
cylinders. Steady state permeation experiments were run over a range of superﬁcial
velocities in both arrays, to study the onset of departure from Darcy’s law. The
rheology of the ﬂuid was evaluated in shear before and after each set of runs.

While departures from Darcy’s law occurred in both arrays at similar values of

102

103
Deborah number, mechanieal degradation of the polymer solution was much more severe

with the triangular pitch array than with the rectangular pitch array. Speciﬁcally, after
several runs through the triangular array the relaxation time was halved while the change
in viscosity was relatively minor; this reveals loss of the high molecular weight tail in
the original polymer. This degradation was irrecoverable; no recovery was noted after
two weeks. Measurements of molecular weight distribution on the same samples in
Odell’s laboratory conﬁrm that the highest molecular weight components are degraded.

Finite element simulations of Stokes ﬂow were carried out for the two different
geometries to determine extensional strain rates along the ﬂow direction in several
regions. This was followed by calculations of polymer chain deformation in these
regions, with the nonlinear elastic dumbbell model. These ealculations reveal that the
maximum stretch rate in the triangular pitch array occurs along the streamline joining the
stagnation points on adjacent cylinders; this leads to nearly complete extension of the
polymer chain at a nominal Deborah number of l in the triangular array. However, in
the rectangular pitch array, the maximum stretch rate occurs along streamlines
considerably removed from the stagnation points, and the polymer chains are not

extended along those streamlines up to a Deborah number of 1.

MW

The ﬂow of liquids through regular arrays of cylinders arises in a variety of

applieations ranging from heat exchangers with tube bundles, to manufacture of ﬁber

104
reinforced composites. Darcy’s law is often employed as a macroscale model for ﬂow

of incompressible, Newtonian ﬂuids through porous media at low Reynolds numbers.
This macroscale representation of the superﬁcial velocity v in anisotropic media can be

written as

Mrs -x-<VP>, (4-1)

where n is the viscosity of the ﬂuid, K is the permeability tensor deﬁned entirely by the
geometry of the array and <VP> is the mean pressure gradient in the ﬂuid. The
longitudinal permeability K33 describes ﬂow along the direction of the aligned axes of the
cylinders, x3. Flow in the plane transverse to the cylinder axis may be described in
general by three constants - K“, K”, Kn, for any given conﬁguration - cf Sangani and
Yao (1988). These constants represent the transverse permeability which is generally
much lower than the longitudinal permeability. In arrays with additional rotational
symmetry such as square or hexagonal packing of cylinders, one parameter sufﬁces to
describe the transverse permeability. Theoretieal values of this quantity have been
tabulated for both square and hexagonal arrays by Sangani and Acrivos (1982). These
values were obtained by numerical solution of the creeping ﬂow equations over
representative cells for these arrays. They have also provided analytical expressions for
this quantity in dilute arrays and in concentrated arrays. For arrays with void fraction
greater than 0.5, the predicted values of transverse permeability with square packing and
hexagonal packing are not signiﬁcantly different.

No quantitative results are available from experiments or from theory for the

105
permeation of viscoelastic liquids transverse to regular arrays of cylinders. The ﬂow of

viscoelastic liquids at low Reynolds numbers through packed beds of spheres has been
studied experimentally by Marshall and Metzner (1967). The ﬂuids they used has a
shearratedependentviscosityasthattheseﬂuidsobeyEq.4.lonlyuptoacertain
critiealvalueofnominalstrainrateinthebed. Asthestrainratewasincreasedabove
this threshold, the pressure gradient or frictional resistance increased progressively from
theDarcyvalueby factors of lOormore. This increasewascorrelated withaDeborah
number, which is the product of a ﬂuid relaxation time and strain rate. Other workers
(James and McLaren, 1975, Kulicke and Haas, 1984 and Durst er al. , 1981) have studied
the resistance to ﬂow of very dilute (”drag reducing”) solutions of polymers passing
through packed beds of spheres. These workers have reported an onset Reynolds number
at which the ﬂow resistance increased suddenly from Newtonian behavior by an order
of magnitude. For example, James and McLaren (1975) worked with dilute aqueous
solutions of polyethylene oxide passing through packed beds at low Reynolds numbers.
They observed that the onset Reynolds number decreased with increasing concentration
and with increasing molecular weight of polymer. They reported also that degradation
occurred especially with larger bead sizes (0.045 cm diameter). Kulicke and Haas (1984)
have shown that, for a given solvent polymer pair, the onset Reynolds numbers may be
used to infer the weight average molecular weight of the polymer. Both effects-the
increase in ﬂow resistance, and degradation-have been attributed to the strong extensional
ﬂow component in such ﬂows (see Durst er al., 1981).

The object of this paper is to report speciﬁc features of this extensional ﬂow

106
component that appear to be critical in crossﬂow of viscoelastic liquids through

unidirectional cylinder arrays with different packing geometries at a given level of
porosity. The work involves both ealculations of chain extension in speciﬁc geometries
and steady state permeation experiments through such arrays. The liquids used have
constant shear viscosity and signiﬁcant elasticity. The results show that the packing
geometry has a signiﬁeant effect on the extent of degradation of the polymer in the
porous medium. This result is explained with model ealculations of polymer chain

extension in the two arrays.

4.3 Experimental

4.3.1 Materials

Three test ﬂuids were used in this investigation-the elastic Ml liquid, another
dilute solution of polyisobutylene in polybutene and a Newtonian analog to the M1 ﬂuid.
It was necessary to prepare another dilute polyisobutylene solution similar to the M1 ﬂuid
because only a limited supply of the M1 liquid was available. This elastic analog was
prepared by ﬁrst dissolving the polyisobutylene (V istanex L-120 from Exxon Chemieal
Company, M. = 1.66x10‘) in kerosene. This solution was then mixed with an
appropriate amount of polybutene (Indopol H300 from Amoco Chemieal Company) in
order to obtain a shear viscosity similar to that of the M1 ﬂuid (see Figure 4.1). The
Newtonian analog was prepared by mixing 17% by weight kerosene into polybutene,
againin order to obtaina shear viscosity similarto thatofthe Ml ﬂuid (seeFigure 4.1).

10

107

 

 

1) [PO'SJ

 

M1 Fluid. 1' - 20 'C

Newtonian Analog. T - 20 'C
Elastic Analog. - 22 'C

 

 

 

2 2 2 2 9 I g
. O 9 O Q . .
T r r r 1 tr . I , r r I I
0.1 1.0 10.0
7" [VS]

Figure 4.1 Steady shear viscosity of the fresh test ﬂuids.

108

4.3.2 Rheological properties

All viscometric measurements were obtained with a Rheometrics RFS-8400 ﬂuids
spectrometer. Viscosity measurements were made under steady shear at rates ranging
from 0.1 to 10 s" and at six temperatures ranging between 15 and 30 °C. The dynamic
moduli of the elastic ﬂuids were measured under oscillatory shear at frequencies ranging
from 0.1 to 100 rad/s. Measurements on the M1 liquid were made at six temperatures,
also ranging between 15 and 30 °C. The elastic analog was tested only at 22 °C.
Figure 4.1 is a plot of the steady shear viscosity, n, of the test ﬂuids. The test ﬂuids
do not exhibit any shear thinning for strain rates less than 1 3". Only slight shear
thinning ean be observed for the M1 ﬂuid and the elastic analog for rates greater than
1 s“. This ﬁgure also shows that these ﬂuids have similar viscosities near room
temperature. This can be seen more clearly in Table 4.1 where the properties of the
ﬂesh test ﬂuids are compared at the same temperature. Figure 4.2 shows the results of
the oscillatory and steady shear experiments, using the M1 liquid, at 20 °C. These data
extend to low enough frequencies and shear rates where the dynamic viscosity matches
the steady shear viscosity. A relaxation time A, may be calculated from the low

frequency region of the storage modulus, G', where G’ is quadratic in to, as follows:

A. =- . .
1 62(ﬂa'ﬂ,) (4 2)

 

Table 4.1 Properties of the fresh test ﬂuids.

 

Fluid

Temp.

 

o 1
'C Pe-s kg/e3 see.
H1 Fluid 22.0 2.73 866. 0.220
Newtonian.Anelog 22.0 2.58 873. ---
838 Polybutane
178 Kerosene
Other PIB Solution 22.0 2.28 871. 0.102

818 Polybutene
18.7. Kerosene
0.25. P18 L-120

 

 

 

 

 

 

110

 

100.007 10
' Temp. - 20.0 -c
- o 6‘ m
o #-
10.00":' D 11' o
I O b
I ’7 o
O
o l-
O

 

1.00 aD.U-u g .0.

Ban
1/ D D D u m
'o D o _
.17

6' [Po]

 

 

0.01 I 1 I T—tTftll T v t I v t I 81" 1

°~‘ 1.0 ' "13.0 100.0
0. ‘i [rod/s, 1/s]

Figure 4.2 Dynamic moduli and steady shear viscosity of the fresh Ml ﬂuid.

[9. 0d] (1. an

111

Here n, is the zero shear viscosity of the solution and n, is the solvent viscosity. This
relations provides a good estimate for the FENE dumbbell model in the form used by
Chilcott and Rallison (1988) if the ratio of fully extended dumbbell length to equilibrium
length is 10 or more for the polyisobutylene is polybutene (cf. Chmielewski at al. , 1990).
Acomparisonoftheproperties fortheMl ﬂuidandtheelasticanalogisalsogivenin

Table 4.1.

4.3.3 Apparatus

A diagram of the experimental apparatus is shown in Figure 4.3. The apparatus
consists of a reservoir connected at one end to a nitrogen cylinder and at the other end
to the permeability cell. During an experimental run, the reservoir is charged with the
test ﬂuid. A constant pressure is then supplied to the reservoir via a nitrogen tank and
pressure regulator. The ﬂowrate was varied by varying the upstream reservoir pressure.
As the ﬂuid exits the reservoir, it passes through a 5.3-to—1 contraction; the ﬂuid then
ﬂows through a 20 cm long pipe to enter the permeability cell with a 2-to-1 expansion.
The permeability cell (with cross-sectional view shown in Figure 4.3) consists of three
sections. The cross-section of the cell is rectangular, measuring 5.33 cm by 1.9 cm.
The enhance section leading to the array is 7.62 cm long. The length of the cylinder bed
is 2.24 cm. This consists of ﬁve rows of cylinders with a maximum of ten cylinders per
row; the cylinders are aligned along the 1.9 cm gap. The exit section is 2.54 cm long.
The test ﬂuids exit the permeability cell by going through a 2—to—1 contraction and into

a short exit pipe.

112

 

 

 

 

 

 

 

 

 

 

n
‘ Reservoir
Pressure
Transducers
Permeability
Cell
(A)

Pressure Taps

 

 

\ _/\___/\_qu

7.88 8.24 8.84

Cell Width: 5.28

(8)

Figure 4.3 (A) Schematic of experimental apparatus and (B) a cross-sectional view of
the permeability cell (all dimensions are in centimeters).

113
Two different cylinder arrays were used in this study as shown in Figure 4.4.

In one array, the rows are aligned one behind the other to produce a rectangular pitch;
in the other array, successive rows are staggered to produce an equilateral triangle pitch.
The cylinder radius a is 0.159 cm and the porosity is 0.704 in both arrays. The gap
between two cylinder surfaces in a row (i.e. perpendicular to the ﬂow direction) is 0.238
cm for either array. The spacing between cylinder axes in successive rows of the
staggered array forms an equilateral triangle. This array is equivalent to the hexagonal
packing referred to by Sangani and Acrivos (1982). However, in the rectangular array,
the gap between cylinder surfaces in successive rows is 0.163 cm. The cylinder ends are
threaded to ﬁt into the top and bottom plates in corresponding patterns. In the rows
containing ten cylinders, the cylinder on either end is tangent to the channel wall.

The pressure drop across the bed of cylinders was measured by two Omega PX-
610 miniature pressure transducers connected to a strip chart recorder. The bed length
I was the same (2.24 cm) in all experiments. The ﬂow rates were measured
gravimetrically. All experiments were conducted at Reynolds numbers (see Eq. 4.4)
ranging from 0.005 to 0.1. The temperature of the ﬂuid was monitored by a small

thermocouple contacting the ﬂuid at the exit of the permeability cell.

4.3.4 Wall and end effects

Larson and Higdon (1987) have analyzed the ﬂow near the surface in transverse
ﬂow of a Newtonian ﬂuid at low Reynolds numbers through periodic arrays of cylinders.
They concluded that the effect of external velocity is damped out well before the second

114

Flow

 

 

 

 

 

 

 

Figure 4.4 Cylinder array geometries: (A) rectangular pitch array; (B) triangular pitch
array (radius = 0.159 cm and porosity = 0.704).

115
row of the array for a solids concentration of 0.3. A crude estimate of the distance into

the array for effective damping (c1: Larson and Higdon, 1986) is given by Vic, where k
is a relevant permeability. Using the transverse permeability values reported in the next
section, we may then estimate that end effects are damped within 0.05 cm. The effect
ofthesidewallsisexpectedtobesmallbeeausetherearetencylindersineveryrow.
Finally, good agreement between the experimentally determined permeability for the
triangle array and the theoretieally value obtained by Sangani and Acrivos (1982)

conﬁrms that the wall and end effects in these experiments are acceptably small.

4.4 W

4.4.1 Transverse permeability of Newtonian liquid

The magnitude of the mean pressure gradient is plotted against superﬁcial velocity
in Figure 4.5 for the Newtonian ﬂuid ﬂowing transverse to the two arrays. The ﬂow
in either bed is seen to follow Darcy’s law. Both lines drawn through the data on this
log-log plot are of slope l and the nansverse permeability may be determined from the
intercepts of these lines. The value of 2.71 x 10" m2 for the array with triangular pitch
compares very well with the theoretieally predicted value of 2.86 x 10” m". The
theoretical value of transverse permeability for the triangular array may be obtained from
the results of Sangani and Acrivos (1982) [in Table 2 of their paper]. The transverse
permeability for the rectangular array with the same porosity is found to be 4 x 10'7 m’.

This is higher than the theoretical value of 2.67 x 10” m2 for a square array with the

116

1 000000

 

 

I Triangular Pitch
D Rectangular Pitch,

 

 

l L

-AP/l [Pa/m]

 

 

 

 

0.01

'o.1o

Figure 4.5 Mean pressure gradient vs. superﬁcial velocity for the Newtonian ﬂuid.

117
same area void fraction presumably because the gap between cylinder surfaces in one

direction (0.238 cm) is considerably higher than the gap in another direction (0.163 cm);

either gap would be 0.2 cm for the square array of the same porosity.

4.4.2 Departure from Darcy’s law for elastic liquids
The onset of viscoelastic effects may be conveniently represented on a plot of the
porous bed friction factor against a Reynolds number. The porous bed friction factor

is deﬁned conventionally by

f: ('AP)P’28€3
162(1-6)

 

(4.3)

where p is the ﬂuid density, -AP is the pressure drop, e is the void fraction of the bed

and G. is the superﬁcial mass ﬂux. The Reynolds number is deﬁned by

am,

Re .3 n,(1-e)'

(4.4)

For Newtonian liquids, the friction factor is inversely proportional to the Reynolds
number. The data for Newtonian liquids ﬂowing through the rectangular array ean be
correlated by the product f-Re = 100 while the data for the triangular array ean be

correlated by an: = 150.

118
The friction factor is plotted against Reynolds number for the M1 ﬂuid through

the triangular and rectangular pitch arrays in Figures 4.6 and 4.7, respectively. Since
only a limited supply of the M1 ﬂuid was available, each of the experimental runs was
made with the same charge ofMl ﬂuid. The ﬁrst set ofexperiments using the fresh Ml
ﬂuidwasmadewiththetriangularpitcharrayandisdenotedasSetlinFigure4.6.
During this set of experiments the superﬁcial velocity varied from 0.17 cm/s to 1.34
cm/s. The onset Reynolds number is 0.01 for the ﬁrst set of runs with fresh M1 liquid.
Asecond setofrunswas madewiththesame ﬂuidinthetriangularpitcharrayisorder
to check for polymer degradation. Figure 4.6 shows that the onset Reynolds number for
this set has increased to 0.02, indicating degradation of the polymer solute. This
observation is consistent with the increase in the onset Reynolds number observed by
James and McLaren (1975) and by Kulicke and Haas (1984) upon lowering the polymer
molecular weight in a polymer solution. Only a few of the data points after the onset are
shown on Figure 4.6 for either set beeause degradation changes the properties of the
ﬂuid progressively during the course of several runs. Data for the rectangular pitch
array in Figure 4.7 were also obtained from runs with degraded M1 liquid and the onset
Reynolds number was around 0.026.

Measurements of relaxation times reported in the next section for the M1 liquid
at various stages of degradation show that the superﬁcial velocity at the onset of
viscoelastic effects seales inversely with relaxation time. This may be restated in terms
of the Deborah number - a dimensionless product of the nominal strain rate (vJa) and

the relaxation time M:

119

 

 

 

 

100000.

. “AAAA‘
3
u a A A
3 \. .
u’ o
c 10000-1 \
O “ o O
2.3 " o o
o It
0: 4
U- 4

A Data 1
rooo~ ° '*"’ “ .. .....T , 1
0.001 0.010 0.1 00

Reynolds Number

Figure 4.6 Friction factor vs. Reynolds number for the M1 ﬂuid ﬂowing through the
triangular pitch array.

120

 

Friction Factor

 

1000 e T T . . . e .
0.01 0.10

 

 

 

Reynolds Number

Figure 4.7 Friction factor vs. Reynolds number for the M1 ﬂuid ﬂowing through the
rectangular pitch array.

121

De 8 ° (4.5)

 

The Deborah number at onset of viscoelastic effects is about 0.3 for the triangular array
and 0.4 for the rectangular array. This may be seen on a plot of the ratio f-Re/(f-Re),.
against the Deborah number in Figure 4.8. The magnitude of the resistance ratio is
consistently lower for the runs through the rectangular array than for the triangular array.
This difference may be due in part to the degradation of the M1 before it is passed
through the rectangular array. Fresh batches of another PIB solution have been used to

explore the extent of degradation alone.

4.4.3 Degradation of elastic liquids

The steady shear viscosity and dynamic shear storage modulus were evaluated for
samples of the M1 ﬂuid taken at three different stages to assess the extent of degradation.
These measurements at 23 °C are reported in Table 4.2. The ﬁrst sample was taken of
the fresh liquid. The second sample was taken after the ﬂuid was passed 20 times
consecutively through the triangular array at nominal strain rates ranging from 0.5 to 8.4
s“. The relaxation time changed from 0.19 s for the fresh ﬂuid to 0.098 s for the second
sample; the shear viscosity change was relatively minor-less than 10% . This trend is
ﬂuid property change is typically associated with changes in molecular weight
distribution. The third sample was taken after the degraded ﬂuid was passed 13 times

122

 

 

 

 

 

 

 

 

10.0. .
3 I Triangular Pitch
~ El Rectan ular Pitch ‘ I a
4 I a
. . . D D D
I“: -l D u
g D
a
53 10 —-——— —D--— D ‘rr 0
\ -l . I—I'I
0 i
Q: .l
J- 1
1
0.1 v . . . . . . .
0.10 ‘ 1.00

Deborah Number

Figure 4.8 Flow resistance of the M1 ﬂuid relative to the Newtonian value vs. Deborah
number.

123

Table 4.2 Properties of elastic liquids before and after crossﬂow runs in different

 

 

 

arrays.
State "a 11
Pa-s sec.
Kl Fluid at 23'C
Fresh 2.53 0.190
After 20 passes through 2.30 0.098
triangular pitch array
After runs through 2.26 0.080
both array types
Two weeks following 2.26 0.080
runs in both array types
Other P18 solution at 22'C
Fresh 2.28 0.102
.After 14 passes through ‘ 2.17 0.065
triangular pitch array only .
After 14 passes through 2.28 0.094
rectangular pitch array only

 

 

 

 

 

124
throughtherectangulararrayatnominalstrainratesranging from 1 to 11 s‘. The

relaxation time changed from 0.098 s for the second sample to 0.08 s for the third
sample with an insigniﬁcant change in shear viscosity. The rheologieal measurements
were repeated after a period of two weeks and the results indieated little recovery.
Molecular weight distributions were determined for the ﬁrst and third samples by Odell’s
groupandarereportedinMiilleretaI. (1990). T'hesemeasurementsconﬁrmthatthe
highest molecular weight components in the fresh sample - with molecular weights of up
to 32 x 10‘ - were lost after these runs. The highest detectable component in the
degraded sample had a molecular weight of 13 x 10‘.

The measurements on Ml reported in the preceding paragraph do not compare
the potential for degradation of identieal fresh batches in either array. This comparison
was made with fresh batches of another polyisobutylene solution. The fresh
polyisobutylene here had a weight average molecular weight of 1.66 x 10‘, leading to a
relaxation time of 0.102 s for the solution. It should also be pointed out that a number
of passes were required to degrade a substantial fraction of the liquid. A fresh batch of
this polyisobutylene solution was passed 14 times through either array at a ﬁxed reservoir
pressure. In the case of the triangular array, the nominal strain rate was around 78 s.
The relaxation time was reduced by 40% under these conditions. In the case of the
rectangular array, the relaxation time was reduced only by 10%. The viscosity change
was insigniﬁcant in both cases. These experimental results established the much greater

extent of degradation in the triangular array for a range of PIB solutions.

125
4.4.4 Analysisofchainextensioninarrays

It is widely recognized (see Ghoneim, 1985, for example) that the converging
diverging geometry of the pores in porous media causes an extensional ﬂow component
thatmaybeassociatedwiththeassociatedwiththeincreased ﬂowresistance for
viscoelastic liquids. However, the strain required for chain scission is typically generated
near stagnation points, where the ﬂuid transit time is highest. Both the arrays in this
study consist of an equal number of stagnation points. However, the experiments clearly
indieate that variations in packing of cylinders at the same porosity and nominal strain
rates have a profound effect on the stretching of high molecular weight polymer chains,
leading to chain scission. Hence, the extent of polymer chain stretching in the two
arrays, especially along the streamlines joining stagnation points, is explored analytically.

The strategy adopted here is that, close to the onset Deborah number, calculations
of chain extension may be based on the Stokes ﬂowﬁeld. Hence, ﬁnite element
simulations were earricd out for two-dimensional Stokes ﬂow within the two arrays to
determine the velocity component it and the stretching component (du/dx) of the velocity
gradient along segments marked in Figure 4.9. The streamline segment joining two
adjacent stagnation points is chosen in either array. The length of this segment within
the rectangular array is only one-fourth of the corresponding segment within the
triangular array. In addition, a segment of the streamline bisecting the gap between the
cylinders is chosen for the rectangular array. These segments are aligned with the bulk
ﬂow direction 1; they coincide with the lines of symmetry in the mesh and with

streamlines in the ﬂow. The only non-zero component of velocity on these segments is

126

 

(A)

 

 

 

 

 

 

(I)

Figure 4.9 Streamlines computed by ﬁnite element simulation of Stokes ﬂow through
(A) the rectangular pitch array and (B) the triangular pitch array (i and ii
denote selected streamline segments).

127
the velocity u in the bulk ﬂow direction. This is plotted for all segments on Figure 4.10.

The velocity varies from 0 to 3 times the superﬁcial velocity over a distance of two
cylinderradﬁalongthestagnaﬁonsueamlinewiﬂrintheuiangulararray. Thereisﬂow
mversalalongthestagnaﬁonsﬁeamﬁnemtherectangularanaywithmtherlow
velocities. The velocity along the bisecting streamline in the rectangular array is much
higher but the variation is signiﬁcantly lower. Figure 4.11 shows that the strain rate
magnimdeswitlﬁntheuiangruaranayvaryﬁomOtoSﬁmesthenominalsUainrate
while they vary from 0 to 1 times the nominal strain rate within the rectangular array.
The latter variation occurs along the bisecting streamline where the velocity is at a high
level throughout. The strain rates along the stagnation streamline within the rectangular
array are less than 10% of the nominal strain rate. To summarize, the maximum strain
rate within the triangular array occurs in the region of maximum transit time. However,
themaximumsuainmtewithindierectangularamyoccursalongﬂiesueunlme
bisecting the gap, where the transit time is low.

We proceed to calculate the chain extension with a FENE model. The elastic
liquidsusedinthiswork maybedescribedbyaFENEdumbbell modelintheform used
by Chilcott and Rallison (1988). This model includes a parameter L which is the ratio
of fully extended dumbbell length to equilibrium length. The extensional viscosity in
uniaxialextensionofthismodeltendstoaﬁnitelimitthatisproportionaltoL’athigh
strain rates (cf Rallison and Hinch, 1988). For example, the value of L = 10 chosen
here would lead to prediction of a uniaxial extensional viscosity that is 100 times the
shear viscosity. Chilcott and Rallison have shown that ﬂow past an obstacle of a Boger

128

 

 

 

 

5.0
-r-— Triangular Pitch
' - - Rectangular Pitch. i
4.0- Rectangilar Pitch. ii

 

 

 

 

 

0.0 1.0 2:0 ' 3.0 t 4.0

Figure 4.10 Velocity proﬁles along selected streamline segments (566 Figure 43)-

129

 

du/dx t (o/v.)

 

 

— Triangular Pﬁch ‘
- - Rectangular Pitch. i
- - Reetan ular Pitch ii

0.0 1.0 2:0 ' 3.0 ' 4.0

 

 

 

 

 

 

Figure 4.11 Strain rate proﬁles along selected streamline segments (see Figure 4.9).

130
liquidismodeledbyvaluesforthispararneterof3-20. Itappearsthatsimilarvalues

model the degradation behavior observed in this work. .
The ordinary differential equations, which describe the evolution of the two
components R1, and R22 of the normalized mean square end-to-end distance R2 of the
polymerchain, aregiveninSection4.5. ThenormalizationissuchthatR’is l at
equilibrium. These were solved with periodic boundary conditions at the ends of the
segments indicated above. A Tchebychef polynomial ﬁt was used for the velocity along
these segments. The results of these computations are plotted for Deborah numbers of
0.4 and 1 in Figures 4.12 and 4.13, respectively. Each ﬁgure shows the ratio of
stretched length to equilibrium length along all three streamline segments. The curves
for the two segments of the rectangular array coincide in both ﬁgures. It is clear from
these ﬁgures that the chain is signiﬁeantly extended only within the triangular array. The
maximumstretch ratioisS foraDeborahnumberof0.4and9foraDeborahnumber
of 1. There is virtually no extension within the rectangular array along either streamline
at Deborah numbers of 1 or less. This is consistent with our experimental observations

on degradation in ﬂow of elastic liquids transverse to aligned cylinder arrays.

4.4-5 W

The force law for the FENE dumbbell is given by

g ;__L__
f(R) R3 (4.6)

131

 

 

 

 

 

 

 

 

12
-— Triangular Pitch _
q -- Rectangular Pitch.i 0° 0'4
10- Rectangular Pitch. ii L - 10.
a-
O’ 'l
as” 5
\.
a: 4
4-4
2-4
0 ' ' Y ' ' ' ' 1 v v v I . r .
0.0 1,0 2.0 3.0 4.0

Figure 4.12 Computed stretch ratio of a polymer chain alongselected streamline
segments: De = 0.4.

132

 

 

 

 

12
-- Triangular Pitch '
. e _ .
-- Rectangular Pitch. i D 1 O
10- Rectamular Pitch. ii L - 10.

 

 

 

 

 

 

Figure 4.13 Computed stretch ratio of a polymer chain along selected streamline
segments: De = 1.0.

133

where R2 is the normalized mean end-to-end distance squared:

R’-R,’+R.’+R§.

(4.7)
At equilibrium
1
R3” —3-. (4.8)
Inplanar ﬂow,
R33 . 3- (4.9)

The evolution of R,‘I and R22 along selected streamlines, where u is the only nonzero

component, is described by the nondimensional differential equations (see Rallison and
Hinch, 1988)

 

 

de 2 du f 2 1 (4 10)
d-R: a du f a 1 (4 11)
U dx ZR; 3 B; (R; 3 )_.

The boundary conditions are periodic:

134

R1(x=0) = Rdxss), (4.12)

thx-o) = R,(x=s). (4.13)

where 0 and s denote ends of the segmt.

4.5 Conehrsians

The ﬂow of dilute solutions of polyisobutylene in polybutene transverse to
unidirectional arrays of cylinders at low Reynolds numbers is sensitive to packing
geometry with a ﬁxed porosity. These liquids have a constant shear viscosity and
relaxation times ranging from 0.1 to 0.2 s. The onset of elastic effects on the permeation
rateordepaereﬁomDarcy’slawoccursatsimilarvaluesofDeborahnumberinthe
two arrays tested here. The extent of degradation of the polymer above this onset
Deborah number, however, is much higher in the triangular pitch array with staggered
rows than in the rectangular pitch array. Calculations of chain extension with a FENE
model of the form used by Chilcott and Rallison reveal that nearly complete extension
of the polymer chain occurs along streamlines joining the stagnation points in the
triangulararray. Thisisbasedonthefactthatthemaximumstretchratewithinthe
triangular array is attained along the streamline joining the stagnation points. However,

in the rectangular array, the stretch rates along the streamline joining the stagnation

135
points are negligible; the maximum stretch rate occurs along streamlines considerably

removed from the stagnation points so that the total strain is insufﬁcient to extend the

polymer chain.

Chapter 5

VISCOSITY EFFECTS IN THE PRODUCTION OF
COMPOSITE PREPREG BY HOT MELT
INIPREGNATION

5.1 Summary

An experimental and theoretieal investigation of the hot melt impregnation process
is discussed. Experimental data show that when a ﬁber tow is pulled through a resin
bathand thenthroughawedgeshapeddie, thetotalresinmassfractionaswellasthe
extent of resin impregnation in the tow depend on the processing viscosity. The
penetration of resin into a ﬁber bundle is greater when the resin viscosity is higher.
Thiseannotbeexplainedbyeapillarity,andthesefeaturesarenotaffectedbytow speeds
up to 25 cm/s. A theoretical model is developed to describe the dependence of
impregnation on viscosity. This model incorporates tow consolidation through both
transverse compression and tension mechanisms. Good agreement with experimental
observation is obtained with a viscosity dependent effective tensile modulus of the

partially wetted ﬁber bundle.

136

137
5.2 W

A key step in the manufacture of complex composite structures is the
impregnation of a bundle of ﬁbers with a liquid resin. The goal of this processing step
is to fully wet all of the individual ﬁlaments which comprise the ﬁber tow, producing
”prepreg” with a predetermined and uniform resin mass distribution. Hot melt
impregnation (HMI) is one such process (Lee et al., 1986). During hot melt
impregnation, a ﬁber tow is unwound from a feed spool, pulled through a resin bath and
then wound up on a large drum (see Figure 5.1 (a)). The actual impregnation of the
ﬁber bundle occurs in the resin pot (see Figure 5.1 (b)) where the tow is pulled through
the liquid resin, past stationary “impregnation” bars and then out through a die located
at the bottom.

The stationary bars have an important role in the impregnation process. Bascom
and Romans (1968) pulled a bundle of glass ﬁbers past stationary bars submerged in
epoxy resin. They observed that as the ﬁber bundle wound around the stationary bars,
resin was squeezed into the tow while air was simultaneously squeezed out, thereby
reducing the tow’s void content. In a study of the HMI process, Chmielewski er al.
(1988) found that when impregnation bars are employed, the resin mass ﬁaction of the
prepreg product is independent of both the processing viscosity and speed. Lee er al.
(1988a) have made similar observations with an identieal HMI process. These results
are in contrast to those obtained when the impregnation bars are not employed. When

prepregging without impregnation bars, lower prepreg resin mass fractions are obtained,

138

 

 

 

I
Resin
Pot [
Wind-up .
Dmm

 

Fiber Spool

 

(a)

l lmpregnation Bars

 

 

 

 

 

( , Exit Die

Fiber Tow

(b)

Figure 5.1 Hot melt impregnation (a) schematic of the hot melt prepregger (b)

schematic of the resin pot.

139
andtheresinmassfmcﬁonisdependurtmﬂwresinviscosityduﬁngpmcessing.
Greater resin viscosities produced prepreg with higher resin mass fractions. However,
the effect of tow speed on the resin mass fraction remains insigniﬁeant (see Figure 5.2).

TheimpregnationbarshaveonefurthereffectontheHMIprocessandthatisto
reduce the maximum tow speed at which the process can operate. Chmielewski er al.
(1988) found that without the bars, speeds of 25 cm/s could be attained, but with the
bars, line speeds greater than 10 cmls result in the ﬁber tow fraying and incipient
breakage at the die. Also, a combination of higher resin viscosities with the bars eauscd
the tow to break at even lower speeds. Using impregnation bars, Lee et al. (1988a) also
found increased tow damage and breakage as larger resin viscosities and line speeds were
employed. They attribute this to increased shear stress and friction within the exit die
alone. However, the above results suggest that greater drag around the impregnation
bars, resulting from higher viscosities and tow speeds, promotes fraying of the ﬁber tow
and then eventual breakage in the die.

Therefore, at fast processing rates that are typieally necessary in order to produce
relatively low cost materials, the HMI process is most efﬁciently operated without the
impregnation bars. This investigation has been undertalmn with the ultimate goal of
increasing theHMIprocessingratebyremovingtheimpregnationbars, andyetstill
producing quality prepreg. This can be accomplished only by understanding both the
impact of processing conditions and the role of the exit die on the impregnation of a ﬁber
bundle.

Experimental work has been done to reveal the distribution of resin within a ﬁber

140

 

 

 

 

 

 

 

0.40
‘ '/ WWW Bor- [Possum-roared]
c I I I g
:2 0.30 I e 0 0 0
O ‘ ° ° e
O . a 0 ° 0 a a
h a . n I g n
L; '1 . D V ' g
to 0.20% a v ‘ ' ‘ ‘ ‘
O a
2 a
c ‘ a a
'63 j
& 0.10-l ' WMPMQ - 0,1 Poo.
. . DGEBA/rnPDA. 1" - 4,5 Poo.
t A 9017501000. ﬂ - 0.8 Pa-s
, U Polybutene n - 4.1 Pa-a
0.00 ' - . . ' . ; . , .
0 5 10 15 20 25

Velocity [cm/s]

Figure 5.2 Resin mass fraction dependence on tow velocity and resin viscosity.

141
bundle when a tow is passed through epoxy resins of widely different viscosities. A

lubrication model is developed to explore the role of the exit die in ﬁber impregnation.
This model incorporates tow consolidation by means of a compressive force due to
hydrodynamic loading and a tensile force resulting from drag on the tow as it travels
through the die.

53W

The effect of processing conditions on the resin content of composite prepreg is
studied via a bench scale hot melt prepregger (Model 30, Research Tool Corporation).
ResinwntentmeasummamaremadeonprepregproceswdmspeedsmngingﬁomOJ
to 25 cm/s and with resin viscosities ranging from 0.1 to 26.6 Pa-s. Resin viscosities
arecontrolledbyregulatingtheresinpottemperamre. Prepregresinmassfractionis
calculated from gravimetric measurements made on tow pieces, 20 cm long, cut from the
tapeexitingtheresinpot. TheresinsusedinthisstudyareaNewtonian,highviscosity
polybutene liquid (Amoco Indopol H300) and an epoxy resin system. The epoxy is
diglycidyl ether of bisphenol-A (DGEBA) and it is mixed in stoichiometric proportions
with m-phenylenediamine (anDA). A 12k tow of Hercules AS4 carbon ﬁbers constitutes
the reinforcement phase.

In order to assess how the resin is distributed, between ﬁber bundle coating and
bundle penetration, photographs of the cross section of the ﬁber tows impregmrted with

the epoxy resin were tam. Two ﬁber tow samples processed at 18.3 cm/s but at

142
different viscosities, 4.5 Pa-s and 0.1Pa-s, were cured at room temperature for 48

hours while being held at the same tension observed during processing. Cross sections
of the cured tow samples were cut and mounted in a quick setting epoxy. These samples
were then polished and photographed at a magniﬁcation of 12.5x.

5.4 Rainmfrastinnandjistrihutiun

The experimental results of Figure 5.2 indicate that the resin content of a prepreg
tape is nearly independent of the rate at which it is processed. However, these same
results show a signiﬁcant viscosity effect on the resin content. As the processing
viscosity is varied from 0.1Pa-s to 26.6Pa-s tbeprepreg resin mass fraction increases
from 19 to 28 percent. Since the data presented in Figure 5.2 offer no information on
the distribution of resin between the surface and the interstices of the impregnated ﬁber
tow, photographs ofthe cross section ofthe epoxy impregnated ﬁber tow were taken (see
Figure5.3). InFigureS.3, thedarkareaswithinthetowsareregions voidofresin
while the white portions of the tow are resin impregnated areas; at this magniﬁcation
individual ﬁbers cannot be resolved. The dark regions surrounding the ﬁber tows,
however, are not void of resin and yet still prevent the measurement of resin coating
thiclmess. However, a comparison of resin penetration can be made. A qualitative
comparison between Figures 5.3 (a) and 5.3 (b) reveals that higher resin viscosities lead
to greater penetration of resin into the ﬁber bundle! The void content of the ﬁber tape

manufactured witharesinviscodty of0.l Pa-sis much higher (darkareasinFigure5.3

143

 

Figure 5.3 Two 12kAS4 carbon ﬁbutowsimpregnated with aDGEBA/mPDA epoxy
resin system under the following processing conditions: tow velocity - 18.3 cm/s (a)
resinviscosity = 4.5 Pa-s (b)resinviscosity = 0.1Pa-s.

144
(b)) than that manufactured using a resin viscosity of 4.5 Pa-s (Figure 5.3 (a)). even

though the capillary number of the former situation is much lower, 0.4, than the latter,

18.0. Thus, eapillarity certainly cannot account for the differences in resin penetration.

5.5 Mathmtisalnadelins

To assist in the understanding of the observed viscosity effect on prepreg quality
a mathematical model is developed, describing the ﬂow processes within the exit die.
A lubrication model is used to describe the liquid ﬂow within the die proper while
Darcy’s law is employed to control the penetration of liquid into the ﬁber tow. Because
the ﬁber tow is not rigid but compliant, its response to external forces is also considered,
and modeled as that of a non-linear spring. The elastic response of a ﬁber network
saturated with resin is an important physical phenomena in all types of composite
processes, and is used in this study to explain how larger processing viscosities result in

prepreg with greater resin penetration.

5.5.1 Tow consolidation models

Fiber tow consolidation is accounted for through two mechanisms: (1) transverse
compression, and (2) tension. The hydrodynamic pressure developed in the die as a
result of its converging geometry will place a compressive load on the ﬁber bundle,
resulting in higher ﬁber volume fractions. Likewise, the pulling force needed to move

the tow through the die at a constant velocity will place a tensile stress on the ﬁbers

145
causing them to straighten and consolidate, also increasing the ﬁber volume ﬁaction.

Gutowski er al. (1987) have modeled compressive loadings on composite
laminates by making an analogy between the deformation of a ﬁber in a network of ﬁbers
and beam deﬂection. This resulted in the following relationship for the elastic

deformation of composite laminates under a compressive load, P.‘(x‘),

l V,(x‘)
V - 1
Pc‘(x‘) = it; °

1 V.
- 1
V,(x‘)

 

‘ a (5-1)

 

ha

where V. is the available ﬁber volume fraction (0.907 for hexagonally packed cylinders),
V, is the inlet ﬁber volume fraction and x.’ is the spring constant. This relationship,
along with the experimentally determined value of x.‘ = 158.6 Pa for well aligned AS4
carbon ﬁbers (Gutowski er al., 1987), is used in this study.

Where the transverse consolidation mechanism provides for non-uniform ﬁber
volume fraction proﬁles, a consolidation mechanism based on ﬁber tension is developed
to provide the inlet ﬁber volume fraction. To develop a tension consolidation model, we
assume that under a tensile force, each ﬁber comprising the tow will support the load
equally. Furthermore, we assume that as the tensile force increases, the tow will
consolidate uniformly. An iterative procedure, described in Section 5.5.3, is employed
to relate the tension on the ﬁber tow to the inlet ﬁber volume fraction. A relationship,

similar to Eq. 5.1, is constructed by making an analogy between ﬁber load bearing and

146
the mechanics of column buckling (Gere and Timoshenko, 1984). Figure 5.4 shows

schematically a portion of a ﬁber located within a network of ﬁbers. In the stress free
state the ﬁber is buckled and males several contacts with surrounding ﬁbers. Ifa length
I of this ﬁber is isolated between adjacent contact points and subjected to a tensile force,
FT', the ﬁber’s response is to straighten. This process is modeled as ifthe ﬁber were
abuckled column respondingtoanaxialtensileloadwhere,

FJ-«l'z . (5.2)

Using similar geometric arguments as Gutowski (1985), the relation between the tensile

force acting on the ﬁber tow and the inlet ﬁber volume fraction is,

 

_1__ -3.

o: o Vain Va

FT ‘11.]. - 1? 1 (5'3)
[7; 7;]

where x,‘ is the effective tensile modulus associated with the tensile force and V... is the
ﬁber volume fraction based on the cross sectional area of the die (V... = 0.65 for these
experiments). The term (l/V... - l/V.) in the numerator has been included to preserve
the conservation of ﬁbers within the die when the ﬁber bundle is free of a tensile force.
Figure 5.5 illustrates the relationship, given by Eq. 5.3, between the inlet ﬁber volume

fraction and the ratio of the tow’s tension to its modulus. It shows that in the limit of

147

 

 

(a) (b) (c)

Figure 5.4 An illustration of the straightening of a portion of a bent ﬁber isolated from
a network of ﬁbers; (a) the buckled portion of ﬁber under consideration in an unstressed
ﬁber network; (b) as the ﬁber network is placed under tension, the ﬁber portion of
interest takes up some of the load; (c) the result is that the ﬁber straightens.

148

0.95

 

I

0.90-

q

0.85-

cl

0.80~
0.75-

0.70-i

d

0.65“

I

0.60J
0.01 0.10 1.00 10.00 100.00 1000001000000

FT/xf

 

 

Figure 5.5 Model non-linear tow consolidation behavior as a result of tension on the
ﬁber tow.

149
small forces, V, approaches Vi, while in the limit of large forces, V, approaches V,.

5 .5 .2 Impregnation model

Figure 5.6 deﬁnes the class of impregnation dies studied. The wedge die
geometry is characterized by thedielength, L, thehalfgap, law” and thediehalfangle
0. The speciﬁcshapeofthedieisgiven by thedistancebetween thediewalland the

cerrterline of the die:

L

‘5"2

 

H+1= tan(0)(1-x)+1,

(5.4)

wherethexcoordinateis scaled withLandHisscaledwith 5W2. Theremaining

variables of interest are sealed in the following way:

P' = Pr. r. = awe/aw.) (5.5a,b)
u' = “U. y‘ = yt’AWz) (5.6a.b)
81' = 8i(%Wz) 82' = 82(‘1‘2Wzl (5.7m)
F' = 2Fr,w,L Q‘ = Q(%W,_)U. (5.8a,b)
r; = r.r, r.‘ = 2r..r,w,L (5.9a,b)

Here the asterisks indicate dimensional quantities, and U, is the tow velocity, W, is the

diewidth,Fistheforceneededtopull thetowatspeed U,and Q'istheﬂuidphase

150

 

 

 

u

Tow Surface "c
i [\l. l ‘7’?“
i l 7'_ l E

 

 

Figure 5.6 Wedge die geometry.

151
ﬂowrateperunitwidth. r,isacharacteristic stress and represents thewallshear stress

for pure Couette ﬂow with a uniform gap KW? A nominal value for r, is 237.2 Pa,
based on a typical HMI line speed of 4.6 cm/s and resin viscosity of 0.41 Pa-s. The
ﬁber tow proﬁle is given by 52(2) and the impregnated region by (s, - s,).

A lubrication approximation (see Dean, 1980, for example) is used to simplify
the equation of motion, yielding the following expression for the axial component of the

resin velocity and the shear stress at the tow-resin interface:

 

 

 

 

1
—w
__y-sz _1_2’dpl._2__
u 1 (H+1)-82+ 2 L d'ﬁl‘" 53) (3’ 82)((H+1)+82)
(5.10)
11W:
_ _ 1 .l 2 d? __ (5.11)
nyIy-I. (H+1)—sz + 2 L dxHHH') 82]

These expressions are derived by satisfying the no-slip conditions at the tow surface and
at the die wall. It should be noted that Beavers and Joseph (1967) propose a slip
boundary condition to be used at the interface of a permeable surface. The work of
Taylor (1971) and Richardson (1971), however, indicates that at the ﬁber volume
fractions of interest here (> 50 %) this condition is unnecessary. The ﬁrst two terms
in Eq. 5.10 are due to Couette ﬂow caused by the moving tow, whereas the third
contribution stems from an induced Poiseuille ﬂow due to the wedge geometry. The

force needed to pull the tow through the die is calculated by:

152

1
p- . f "rely-u. dx (5.12)
0

Mass balances on the resin within the die gap and the tow yield the following

results:

roll-

W3 io‘f—up, forsl>o (5.13)
L dx i O. forslso

 

1
3 ”a
L

d _[+u . fors >0 (5.14)
E[(Sz-81)(1‘Vf)] "l g , for 81"0

In the above expressions, s10!) represents the resin impregnau'on front within the ﬁber
bundle, where s,(x) = 0 indicates a fully impregnated tow, and Vfoc) is the ﬁber volume

fraction. The resin phase ﬂow rate is represented by:

3’1

0 - fulx.yldy (5.15)
a,

The normal component of the resin velocity, u,,, at the tow interface is related to the

local transverse pressure drop within the tow by using Darcy’s law:

153

P-( ,-P)
up.xr{ 3:3, } (5.16)

 

In theaboveexpression Kris thetransverse permeability, l‘ is tlreresinpressure dueto
eapillarity and p, is the atmospheric pressure. The longitudinal ﬂow of resin within the
ﬁber bundle as the result of a local axial pressure gradient will be neglected. Only the
case where the tow velocity is much greater than the longitudinal Darcy velocity is
considered.

The Blake-Kozeny-Carman equation (Bird er al. , 1960) is used to relate the

transverse permeability to the ﬁber volume fraction:

(I'V:)3 1 2

V: 4k(V:) I! (5.17)

KTI'

where the dependence of the Kozcny constant, kW), on the ﬁber volume fraction is
based on model calculations of a Newtonian ﬂuid ﬂowing through an inﬁnite array of
circular cylinders (Sangani and Acrivos, 1982). The ﬁber radius is r,

The capillary pressure is given by:

154

P = 2“;:::1’) (5.18)

The resin-ﬁber contact angle is represented by a, and Na, is the capillary number:

- "U° (5.19)
Yr!

 

where his the resin-ﬁber surface tension. Variables found in Eqs. 5.16-5.19 are scaled

as follows:

Kr' = 190151112)2 (51kb)

up. = “DU.

1" = Pr, r,’ = rK‘AWz) (521a,b)

For these calculations 7,, is taken to be the surface tension of DGEBA (44 dynes/cm at
60 °C) as determined by the pendant drop method (Weaver, 1982). It is more difﬁth
to ﬁnd a value for ¢, however. Because impregnation ﬂows are non-equilibrium
prosesses, the dynamic contact angle, and not the static one, is more important when
describing the advancing liquid ﬂow front. Inverarity (1969) studied the wetting of glass
and polymeric ﬁbers with viscous liquids and found that the advancing contact angle is
a complex function of both the ﬁber velocity and the ﬂuid viscosity. We will take 35°
to be the nominal value of ¢ (Lee et al., 1988b), and note, that although Ahn er al.

155
(1991) found signiﬁcant interfacial effects when impregnating a woven fabric

reinforcemert at high ﬁber loadings (V, > 0.50), the capillary pressure had little
inﬂuenceon the results of these model calculations at the NO, levels studied here.
Eqs. 5.13-5.19 constitute the physieal model used to ﬁnd the pressure,

impregnation and ﬂow rate proﬁles within the die. To this end, Eqs. 5.13-5.15 are

 

 

reanangedas:
732;“: 11' B(x)[20- [(H+1)‘52]} (5.22)
5”:
d _ L . - .. _
33-2 1WK,.(1 Vt) [P (P. Pl] (5,23)
- 2
2
Q=Qe+ﬂ;-‘/3 (5.24)
where
B”) ‘5 (5.25)

 

[(31’1) '52]3

In Eq. 5.23 the effective impregnation is represented for convenience by the quantity:

156
s= [(83-81)(1-V,)]’ (5.26)

Eqs. 5.22-5.24 are analyzed according to the following boundary conditions:

(i) P=P,=0atx=0 (5.27)
(ii)S=S,=0atx=0 (5.28)
and
(iii)P =p. =0 at x =1 (5.29)

Neglecting the liquid head in the resin pot, no net pressure drop will be imposed across
the die. Also, for the situations studied experimentally in this paper, Na, > 0.4 and the
tow residence time in the liquid bath is at most three seconds. Hence, it will be assumed
the tow entering the die contains no resin. This boundary condition allows for the study

of the factors which affect the impregnation process within the die proper.

5.5.3 Solution procedure
These equations are solved by a semi-implicit Runga-Kutta method with variable
step size (hiichelsen, 1976 and Villadsen and Michelsen, 1978). At x = 0, the
conditions for P and S are provided, and an initial guess for the inlet ﬂow rate, 9,, is
obtained from the non-porous situation (Middleman, 1977). Also, initially a constant tow
proﬁle is assumed, s2 = s2, for all x. A shooting strategy is employed to reach the
boundary condition at x = l. A means of updating Q, for the ensuing iteration is

obtained by integrating Eq. 5.22 over the die length and then rearranging the result.

157
Msleadstomefouowingintegralpropertyfora:

 

0. = , 1 [ism (JS-Jggldx + éiem [(H+1>-s,1dx]
IBM) dx ° o (5.30)
O

The calculation strategy entails three nested iteration loops in order to ﬁnd the
pressure, impregnation, ﬂow rate and tow proﬁles. During the inner iteration loop, a
constant tow proﬁle s,(x) is assumed and Eqs. 5.22-5.24 are repeatedly solved to ﬁnd
the correct inlet ﬂow rate. The pressure proﬁle resulting from this calculation is then
used to evaluate the transverse consolidation. Eq. 5.1 utilizes the loeal hydrodynamic
pressure in the die to calculate the loeal ﬁber volume fraction. The revised tow proﬁle,

s,(x), is found by recognizing that the cross sectional area of the tow ﬁbers are

conserved:

V0
5,00 = 3,00 + W (520’51(x” (5.31)

Subsequently, this updated tow proﬁle is used to ﬁnd new pressure, impregnation and
' ﬂow rate proﬁles.

In the outermost iteration loop, the ﬁber tension is determined by Eqs. 5.11 and

5.12, and Eq. 5.3 is used to amend the inlet ﬁber volume fraction, V,. Then a revised

158
inlettow positionisfoundbyagainrecognizingthatthccross sectionalareaoftheﬁbers
passing through the die is conserved:

 

(5.32)

5L6 gkmnuaﬁenrenuhuuuthrmsmni

Figures 5.7 (a), (b) and (c) show several pressure, impregnation and flow rate

proﬁles based on the following dimensionless parameters:

a = 15° umwg = 69.2 (5336.6)
r/(‘AWQ = 0.05 v, = 0.76 (53m)
v, = 0.90 x. = 0.67 (5356.6)

Theseparametersconespondtothedesignandoperaﬁngcondiﬁonsusedinﬂte
experiments of Figures 5.2 and 5.3. (Although the correct shape of the experimental
exit die is a wedge followed by a short rectangular slit region, previous calculations
indicate that the rectangular region has little effect on the computed profiles.) The actual
die dimensions and ﬁber radius are:

159

 

 

 

 

 

 

 

 

 

 

80
d — ‘v - 0.001
70‘ --- ‘1’ - 0.010 .
‘ ---- x, - 0.100 °-
...... 3' - 1.000 ...
60"---«:,-10.000 -'
50—
o. 40- 3"
30- 5‘.
.0
. :4
a, '
20" '1
10"
q
O‘WA I J¢ ; # l A ‘ 1 r
0.0 0.2 0.4 0.5 0.8 1 .0

Figure 5.7 (a) Calculated pressure proﬁles for several different values of the effective
tensile modulus.

160
W, = 0.0159 cm L = 0.550 cm (53am

W, = 0.5588 cm r, = 3.9 x 10" cm (S.37a,b)

Figure 5 .7 (a) shows that the peak pressure is strongly dependent on the ﬁber
bundle’s effective tensile modulus. As it, is increased from 0.001 to 10.0 the peak
pressure increases from 20 to 70 times the characteristic stress. Over this range of x,-
the ealculated dimensionless tensile force on the ﬁber tow changes from 0.60 to 0.87 (8.7
to 12.7 mN or 2.2 to 3.2 psi for a 12k ﬁber tow moving at 4.6 cm/s through a 0.41
Pa- 3 viscosity resin). Thus, as x,- is increased by 4 orders of magnitude the ratio F1/x,
decreases nearly the same amount and, as is shown in Figure 5.5, the inlet ﬁber volume
fraction decreases. This leads to smaller clearances between the tow and the die (Eq.
5.32). And, similar to the non-porous lubrication ﬂows (Middleman, 1977), the peak
pressure in the die is sensitive to this gap, where smaller gaps result in larger pressures.
The pressure proﬁles of Figure 5.7 (s) also show that the greatest contribution
of thepressure loading on the tow occurs in the last 15 percent of the die. As a result,
tow consolidation due to the transverse loading has an insigniﬁcant effect on the
calculated resin mass fraction of the tow. For example, for the situation where the peak
pressure is the greatest, x,- = 10.0, the outlet ﬁber volume fraction is less than 8 percent

greater than the inlet value.

0n the other hand, the impregnation proﬁles of Figure 5.7 (b) show that the

161

 

 

 

 

 

 

 

 

 

1.0
—— z, - 0.001 .
4 --- s, - 0.010 . .'
---- e, - 0.100 . - .-
g 0-3" ------ x, - 1.000 . .
:5 ' ' ° r. - 10.000 , . ' .......
o 4 , - .....
c o . . .......
8‘ 0.64 .- .......
L. . o ....
o. ‘ - ..-—‘
E - ' . ........
a) 0.4- . ' , ————
.2 ' .... "’ a
u . ', .—r
o .1 a... ”’o
3 I,“ .—
05 0.2- x” ____________
00. ” .... ﬂﬂﬂﬂﬂﬂ .-
1"; ”M-
I
0-0 I I I ' I I r I '
0.0 0.2 0.4 O 6 O 8 1 O
X

Figure 5.7 (b) Calculated impregnation proﬁles for several different values of the
effective tensile modulus.

 

162
calculated ﬁber tow penetration is sensitive to ﬁber tension effects. The increase in tow

impregnation with x, at each location along the die can be explained by again noting the
decrease in V, withincreasesinx, Astheﬁbervolumefractiondecreases,boththe
transverse permeability of the tow and the hydrodynamic pressure in the die will
increase. According to Darcy’s law (Eq. 5.16), this will result in higher resin penetration
rates.

Figure 5.7 (c) presents the dimensionless resin phase ﬂow rate proﬁles. At each
location along the die the ﬂow rate decreases with increasing x1, opposite to the trend in
FigureSJ (b), inordertopreserve themassbalance. ThedecreaseinQalongthe
length of the die for each value of It, reﬂects the loss of ﬂuid from the resin phase to the
tow. However, atx = 0.97 for the x,- = 10.0 curve, the ﬂow rate proﬁle changes from
a decreasing function of x to an increasing function. At this location, the tow is fully
impregnated and yet still is being consolidated by the hydrodynamic pressure. The tow
consolidation causes ﬂuid to move from the ﬁber bundle into the resin phase, resulting

in an increase in Q.
Figure 5.8 (a) shows model calculation results on how the resin coating and resin
impregnation depend on the effective tensile modulus. These results are combined in
Figure 5.8 (b) which presents the predicted resin mass fraction versus the effective
tensile modulus. Varying the modulus from 0.06 to 10 results in an increase in the
predicted resin mass fractions from approximately 20 to 27 percent. This is similar to
the mass fraction changes observed experimentally by varying the resin viscosity from

0.1Paos to 26.6Pa- s (Figure 5.2). Moreover, Figure 5.8 (a) shows that this increase

163

 

 

 

 

 

 

 

0.4
. — x, - 0.001
b --- It, - 0.010
. ---- at, - 0.100
"3:. ''''' I, - 1.000
0'3 . - , e, - 10.000
§~1e___
‘ sg‘ \.' ____________________
«I ‘~~“~‘h~ .....
or 0.2- ' ""31 -----------
# . O ................
0.1- ' ' - . . .........
O-O ' I I I r I ﬂ r '
0.0 0.2 0.4 0.6 0.8 1.0

 

Figure 5.7 (c) Calculated ﬂow rate profiles for several different values of the effective
tensile modulus.

'164

 

 

 

 

0.30‘ FLO
0. 5'1
2 “0.8
at
g 0.20-
3 r -0.6
O
o o.15-
C
'6 -o.4
6‘2 0.10-
0.05-' h“0.2
Oomlm_0.0
10" 10" 10" 10" 10" 10 10'
E

Figure 5.8 (a) Calculated resin distribution in a mode] ﬁber tow.

uogoubeder entropy

 

Resin Mass Fraction

165

 

0.30
028-1

026-1

.4

0.24-

q

0.224

J

0.20-j

10"

 

10"

10"

10"

10"

10'

10‘

 

10'

Figure 5.8 (b) Calculated resin mass fraction of a model ﬁber tow. -

 

166
is a result of greater impregnation of the ﬁber tow. Again, this same effect is seen

experimentally when the viscosity is increased (Figure 5.3).
The simplest model relating the effective ﬁber tow tensile modulus to the viscosity

canbeconstructedbyamixingrule:

‘1' 3 V:.Gp + (1'Vg) .G1 (5.38)

where G, and G, are the partially wetted ﬁber tensile modulus and resin tensile modulus
respectively. The resin modulus should be approximately proportional to the resin
viscosity and, based on the experimental observations, the resin volume fraction will be
an increasing function of the viscosity. The experimental results suggest that the
effective tensile modulus of a partially wetted ﬁber tow varies with viscosity in the

manner shown in Figure 5.9.

5.7 (landmine

The effect of processing conditions on the resin mass fraction of prepreg tape
produced by hot melt impregnation without the assistance of impregnation bars is studied.
Experiments performed with polybutene and epoxy resin show, not only an increase in
the prepreg resin content with increasing resin viscosity, but also greater penetration of
resin into the ﬁber tow with larger processing viscosities. These observations cannot be

explained by capillarity. A lubrication model of a wedge shaped impregnation die is

167

 

 

 

 

 

10.00: -
i
.1

1.00-i

:2" 1

1

0,103
a

0.01 1 . ...T... . 1 . ..T..
0.1 1.0 10.0

Viscosity [Po ~s]

_ Figure 5.9 Relationship between the partially impregnated resin tow tensile modulus and
the resin viscosity.

168
developed which incorporates ﬁber tow consolidation through a ﬁber tension

mechanism.This model is able to predict the experimental observations with a tensile

modulus that increases with resin viscosity, in a partially impregnated ﬁber tow.

Chapter 6

CONCLUSIONS AND RECOMNIENDATIONS

6.1 mm

The ﬂow of viscous and viscoelastic liquids transverse to periodic arrays of
circular cylinders is studied. Two array geometries, square and hexagonal, are
examined; each having a void fraction of 70 percent. The elastic ﬂuids include several
non-shear thinning dilute polymer solutions composed of a wide range of molecular
weight polyisobutylenes (0.9 to 6 million) dissolved at a concentration of 0.20 wt. % in
a purely viscous polybutene solvent. These ﬂuids are standard type Boger ﬂuids. In
addition, a shear thinning 2 wt. 96 polyisobutylene in decalin solution is studied.

At a void fraction of 70 percent Stokes ﬂow simulations indicate little difference
in the ﬂow resistance between periodic square and hexagonal arrays. Using purely
viscous liquids, these predictions are conﬁrmed here for Reynolds numbers less than
0.30.

In this study particular attention is devoted to the inﬂuence macromolecular
conformation on the enhanced pressure drop across the arrays. The Boger liquids used

in this study are formulated to produce solutions with different levels of relative

168

 

 

..n

 

169
macromolecular extensibility. Being theta systems, the square of the degree of polymer

extensibility associated with each of these solutions is proportional to their molecular
weight. This is consistent with extensional viscosity measurements made by ﬁber
spinning.

The initial departure from Darcy’s law for the non-shear thinning solutions is due
to enhancement in ﬂow resistance and occurs at Deborah numbers of 0.80 and 0.35 for
the square and hexagonal arrays respectively. These onsa values are independent of
molecular weight. Thus, the inﬂuence of molecular weight on the onset of elastic effects
can only be accounted for through the solution’s relaxation time. The differences in the
onsetDeborahnumbers between thetwo array typesaretheresult ofdifferentdegrees
of stretch rates occurring at a constant superﬁcial velocity. By redeﬁning the Deborah
number in terms of the maximum anay extension rate, as determined by the Stokes ﬂow
ﬁeld, instead of the average strain rate, the onset Deborah number in both arrays collapse
to approximately 1. This Deborah number corresponds to the coil-to-stretch transition,
whereamacromoleculeundergoesanabrupttransitionfromacoilcdstatetoan
elongated uncoiled state. This provides evidence that extensional ﬂow effects determine
the onset of elastic effects and dominate the ﬂow processes at high Deborah numbers (De
> 1).

At a constant Deborah number above the onset, the ﬂow resistance is greater for
higher molecular weights, and is consistently higher for the hexagonal array than the
square array. At large values of Deborah number the relative ﬂow resistances in both

arrays become independent of Deborah number. The magnitude of these asymptotic

170
values are proportional to the molecular weight. Differences in the relative ﬂow

resistance at De > 1 when plotted versus the Deborah number indicate that the solution’s
relaxation time alone is not sufﬁcient for parameterizing the enhanced pressure drop.
The inﬂuence of the relative polymer extensibilities of each solution must also be taken
into account. The asymptotic ﬂow resistance - molecular weight relationship is consistent
with the extensional viscosity measurements since the transverse cylinder array ﬂows are
dominated by extensional ﬂow effects at high Deborah numbers. Thus, the experimental
results follow the expected trend of the asymptotic ﬂow resistance with molecular weight
for theta systems.

Measurements of the extent of polymer degradation for non-shear thinning elastic
ﬂuids in rectangular and hexagonal arrays are also studied. The extent of degradation
in both arrays is characterized by large reductions in the polymer relaxation time after
the ﬂuid passes several times through the arrays. At a comparable number of tests, the
amount of degradation observed in the hexagonal array is much more severe than in the
rectangular array. Further evidence that molecular extension in the hexagonal array is
greaterthanthatinﬂierectangularanayisgivenbychainextensioncalculaﬁonsbascd
on the Stokes ﬂow ﬁeld.

The initial departure from Darcy’s law for a shear thinning solution passing
through a cylinder array is a reduction in the ﬂow resistance. This is eventually followed
by an enhancement as the Deborah number is increased, and occurs at Deborah numbers
of 4.2 and 3.0 for the square and hexagonal array respectively. These data illustrate

another important feature of ﬂuid rheology on the transverse ﬂow past cylinder arrays.

171
The ﬂow dynamics are dominated by shear effects at low to moderate Deborah numbers

which result in a reduced pressure drop across the array for shear thinning ﬂuids. At
high Deborah numbers the extensional effects become important, however.

The ﬂow kinematics in both arrays is elucidated by laser Doppler velocimetry and
streak photography. These techniques reveal a ﬂow transition from steady to unsteady
motion at Deborah numbers corresponding to the onset ofenhanced pressure drop for the
non-shear thinning elastic ﬂuids. These results, along with the failure of current steady
state ﬂow calculations to predict the large experimentally measured ﬂow resistances,
imply that proper mathematical modeling of these type of ﬂows must take into
consideration the viscoelastic transition from steady to unsteady ﬂow.

Finally the effect of shear viscosity on the permeation of a viscous liquid through
a compliant array of cylinders is studied. This is accomplished by an experimental and
theoretical investigation of the hot melt impregnation process. It has been found that by
pullingaﬁbertow throughaviscousresin, andthenthroughawedgeshapeddie, the
amount of material impregnating the tow is a function of the resin viscosity; the higher
the resin viscosity, the better the penetration of resin into the ﬁber tow. A mathematical
model, incorporating two ﬁber bundle consolidation mechanisms, results in good

agreement with experiments with a viscosity dependent partially wetted tow modulus.

172

mm

Thisshrdyisaﬁntstepinundershndingtheﬂowofviscoelasﬁcﬁquidsﬂowing
transverse to arrays of cylinders. Its motivation has been to identify and solve some of
the problems associated with processing composite materials containing a viscoelastic
component. Continuing on toward this goal, the next phase of this study should include
conducting crossﬂow experiments with model arrays at smaller void fractions (most
composite materials have void fractions of 30 to 40 percent rather than 70 percent).
Next, the hot melt impregnation study of Chapter 5 presents interesting and counter-
intuitive results in that higher resin viscosities lead to better impregnated ﬁber tows.
This may be the result of tow consolidation effects caused by the tensile load placed on
theﬁbers. Anexperimentalstudyrelatingtheﬁbertensiontoprepregqualityis
suggested. This should accompany consolidation experiments using model cylinder
arrays and both viscous and viscoelastic liquids.

Finally, the unsteadiness of the viscoelastic ﬂows has not sufﬁciently been
explored. The implication that theoretical study of viscoelastic ﬂuids ﬂowing transverse
to cylinder arrays must involve unsteady state calculations should make further study of
the phenomena a high priority of non-Newtonian ﬂuid machinists. The investigation of
the unsteadiness should begin with the analysis of the ﬂuctuating downstream pressure
signals for high Deborah number ﬂows. Next, the problems encountered with the LDV

system (speciﬁcally the poor signal quality) must be resolved. Suggestions for modifying

173
the LDV set-up to obtain quantitative information about the ﬂow unsteadiness include the

following items. The use of new optics to reduce the size of the measuring volume.
This would reduce the number of scattering particles in the measuring volume and reduce
the fringe spacing, assisting in the measurement of low velocities. Furthermore, the
cylinder arrays and windows should be made of optical glass to reduce unwanted light
scattering. Also, the index of refraction should be matched between the test liquid and
cylinder arrays in order to be able to map the velocity in the entire ﬂow ﬁeld without
being hindered by reﬂections from curved surfaces. Lastly, a frequency analyzer is

needed to obtain the power spectrum of the unsteady ﬂows.

APPENDIX

Appendix A

WEDGE DIE/FIBER IMPREGNATION

FORTRAN PROGRAM

A.l Emma < Stlart j

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(Input Data
Initial 82(0)
' O
1
l
82(x)-82(0)
<2)
Q(O)non-por.
<0

 

63

Figure A.1 Schematic of the computer algorithm.

- 174

175

Integrate to find
P00. $00. 000

 

 

 

 

 

 

01°)...” -f(P,S)

 

 

 

 
     

remnant”: < to

Figure A.l Schematic of computer algorithm (continued).

176

Tow Consolidation
V,(X) " 1‘(F’(X))

 

 

 

 

 

 

82(x) - f(S1(x),Vf(x))

 

 

 

 
    

"32 '00-82” (x)|: < to!

Figure A.l Schematic of computer algorithm (continued).

177

 

Tow Consolidation
V, (O) - f(Force)

 

 

 

 

32(0) ' f(V,(0))

 

 

 

 
   

:sz' (0) - 8251(0): < tol

 

 

yes
C End )

Figure A.I Schematic of computer algorithm (continued).

A.2

10.

ll.

12.

13.

178

Routing

. PREPREG

AJACOB

BETA

BJACOB

CONSLP

CON SLT

CSFIT

DFDX

DWALL

ERROR

FLOW 1

FLOW2

FORCE

auction

Main Program
Calculate the Jacobian of
the equation set

Miscellaneous functions

Modification of the
Jacobian

Consolidation of fiber
tow due to hydrodynamic
pressure

Consolidation of fiber
tow due to ﬁber tension

Cubic spline fit of tow
profile and fiber volume
fraction proﬁle

Derivative
set

of equation

Slope of die wall

Maximum difference
between full and half
steps during integration

Inlet flow bounds for
initial guess of Q(0)

Integral property for the
update of Q(0)

Calculate force needed to

pull tow

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

26.

27.

23.

179

Routine

FUNCT

KOZ

Kt

MATINV

MATMULT

RUNGKUT

F 52

SEVAL

SIGMA

SIMPSON

SPLINE

START

STEP

FV

WALL

Bustier
Equationset
Calculate the Kozeny
constant
Calculate the transverse

permeability by the
Blake-Kozeny-Carman
equation

Invert 2 by 2 matrix

2 b y 2
multiplication

matrix

Semi-implicit
integration

Runga-Kutta

Value of the fiber tow
boundary position at
locationx

Cubis spline intepolation
Non-linear spring law
describing tow during

transverse consolidation

Integation
rule

by Simpson’s

Cubic spline ﬁt

Euler scheme to start

integration
Calculate next step size

Value of the ﬁber
fraction at location x

volume

Calculate the location of

the die wall

180

A.3 ﬂagrant
c prepreg
c updated: 02 May 1991
c craig chmialawski
c sunanary: lubrication model of a Newtonian fluid impregnating
c a continuous fiber tow as it is being pulled through
c an exit dis. Darcy's law models the impregnation and
c the tow is modeled as a non-linear elastic network.
implicit integer (i,j,n)
implicit double precision (A-H,K-M,O-Z)
parameter (inc-50001)
dimension Ptinc), Q(inc), S(inc), s2(inc), Vtinc), x(inc)
dimension arr(2), y(2), y0(2), yl(2), y2(2)
common c, DELTAP, GAMMA, kappaP, kappar, Patm
common Pin, 00, rt, 80, s20, theta, V0, Va, Vt
common/ERROR] Btol
common/CSPL/ xxtinc), 822(inc), VV(inc)
common/CSPLsZ/ Bs2(inc), Cs2(inc), Ds2(inc), ndatsz
common/CSPLVI BV(inc), CVtinc), DV(inc), ndatV
external CONSLP, CONSLT, CSFIT, ERROR, FLOII, PLOWZ
external FORCE, Ps2, EV, RUNGKUT, STARE, STEP
naorax(dsg) - deg*pi/180.
c ***** input dimensionless parameters *****
open (unit-1,£ila-'in.put',status-'old'i
raad(l,*) theta, c
read(l,*) rt, Vt, Va
read(1,*) kappaP, kappa?
raad(l,*) Ca, phi
raadtl,*) Pin, Pout, Patm, pct
read(1,*) rhor, rho!
raadt1,*) Ctol, Etol, Qtol, rtol
close (unit-1)
c ***** mdsc parameters *****
pi - 4.*oarau(1.o+000)
theta - RADIANttheta)
phi - RADIAN(phi)
GAMMA - 2.*pi*DCOS(phi)/(Ca*r£)
DELTA? - Pout - Pin
c ***** initial tow-die clearance: s20 - 2.*rf *****

320 - 1. - 2.*rf
V0 - Vt/s20

10

20

3O

35

181

iterT - 0
iterr - iterr + 1
write(*,999) iterr

s10 - s20*(1. - pct/100.)
SO - DSQRT((sZO-s10)*(1.-V0))

***** initial s2 and V *****

inc0 - 10

do 10, i-1, inc0 + 1
x(i) - (i-1)/DFLOAT(inc0)
s2(i) - s20

V(i) - V0

n - incO + 1

iterc - 0

call CSFIT(n,x,s2,V)
iterc - iterc + 1
write(*,1000) iterC

***** initial flow rate guess *****
call FLON1Cle,Qub)

00 - 01b

iterQ - 1

RES - 1.

write(*,1010)

it (mascara .gt. Qtol) then

***** boundary conditions *****

n - 1

x(n) - 0.D+000
P(n) - Pin
8(n) - $0

0 (n) - 00
s2(n) - s20
V(n) - V0

dz - 1./1000.
***** start integration *****

do 35, i-2, 3
3(1) - x(i-1)
P(i) - P(i-1)
8(1) - 8(1-1)
call START (dx,x(i) ,P‘i) ,S(i) ,Q(i))
82(1) - Fs2(x(i))
V(i) ' FV(X(1))

182

ttttt main *tttt

n - 3
:0 - x(n)
y0(1) - P(n)
y0(2) - S(n)
40 if (:0 .lt. 1.) then

***** init full step calc *****

do 50, i-l, 2
50 y(i) - yOCi)

call RUNGKUT¢X0,dx,y)

do 60, i-l, 2
60 y1(i) - y(i)

***** half step calc *****

qq - 1.
70 it (qq .ge. 1.) then
do 80, i-1, 2
30 Y“) ' YOU.)
:1 - 80
dx - 0.5‘dx

do 100, i-1, 2
call RUNGKUT(xl,dx,y)
it (1 .eq. 1) then
do 90; j-1' 2

so we» - yd)
else
endit

100 :1 - :1 + d:

it (y(l) .eq. 0.D+000 .or. y(2) .eq. 0.D+000) then
err(l) - 0.0D+000
err(2) - 0.0D+000

qq - 2.5D-001
else

call BRROR(y1,y,err,qq)
endif

if (dz .lt. 1./FLOAT(inc-1)) qq - 2.50-001

do 110, i-l, 2
110 y1(i) - y2(i)

goto 70

else
endit

183

***** reinitialize *****

dx - 2.*dx
:0 - x0 + dx
do 120, i-l, 2
120 y0(i) - y(i) + 1./7.*err(i)
if (y0(1) .lt. 0.D+000) y0(1) - 0.D+000

***** adjust step size *****
call STEPtqq,dx)
***** final step: adjust step to die exit *****
if (30+dx .gt. 1.) dx - 1.-xO
***** full impregnation limitation *****
if (DSQRI(S(n)) .ge. s2(n)*(1.-V(n))) then
y0(2) - (Fs2(x0)*(1.-FV(80)))**2
else
endif
n - n + 1
x(n) - 30
Pm) -y0(1)

8(n) - y0(2)
Q(n) - 00 + DSQRT(SO) - DSQRI(S(n))

s2(n) - Ps2(x(n))
V(n) - £1”me
goto 40

else

call PLOW2(n,x,s,RES)
write(*,1020) iterQ, 00, RES

endif

***** convergence (regula falsi) to proper QO *****

if (iterQ .eq. 1) then

A - 00

PA - R38

00 - Qub
elseif (itero .eq. 2) then

B - 00

EB - R88

00 - (A*FB - B*FA)/(FB - EA)
elseif (EA*RES .lt. 0.) then

B - 00

PB - R38

00 - (awe - B*FA)/(FB - FA)

130

140

184

else

a - co

m - R88

00 - (A*FB - B*FA)/(FB - FA)
endif

itero - iterQ + 1
goto 30

else
endif

***** tow consolidation by pressure *****

do 130, i-2, n
if (9(1) .gt. P(i-1)) then
Xnax - x(i)
Pmax - P(i)
call CONSLP(P(i),S(i),s2(i),V(i))
else
s2(i) - s2(i-1)
V(i) - V(i-l)
endif
continue

if (iterc .lt. 10) then
do 140, i-2, n
conerr - DABS((Fs2(x(i))-s2(i))/s2(i))
if (conerr .gt. Ctol) goto 20
continue
else
endif

tiiit c‘lc forc‘ *tttt

print*, ' '

print*, 's20- ', s20

print*, 'VO - " V0

print*, ' '

call FORCB(n,x,P,s,forc,1oad)

**‘** tow consolidation by tension *****

s200 - s20

call CONSLT(forc)

ptiht‘, '820- 'g 320

print*, 'V0 - ', V0

print‘, ' '

if (DABS(s200-s20)ls20 .ge. Ttol) goto S

***** calc coating thickness *****

coat - 0(a)

150

185

***** calc fiber vol fract *****
Vf - s20*V0/(s2(n) + coat)
****‘ calc fiber mass tract *****

denratio - rhor/rhof
Hf - s20*V0/(s20*V0 + (DSQRT(S(n)) + coat)*denratio)

eeaaeattaeeeeeeeeeee OUTPUT aeeeeeeeaeeeteetttee

open(unit-1,fi1e-'out.dat',status-'nev')
openIunit-2,file-'P.dat',status-‘new')
open(unit-3,file-'Q.dat',status-'new')
open(unit-4,file-'impreg.dat',status-'neu')
open(unit-5,file-'s2.dat',status-'nev')
open(unit-6,file-'V.dat',status-'new')

write(1,1030) 'theta (deg) ', theta*180./pi
write(1,1030) 'L/(1/2*W2) ', c
write(1,1030) 'rf/(1/2*W2) ', rf
write(1,1030) 'Vf init ', Vt

write(1,1030) 'kappaP*(1/2*l2)/(eta*00)', kappa?
write(1,1030) 'kappar*(1/2*l2)/(eta*00)', kappa!

write(1,1030) 'Ca - (eta*00)/gama ', Ca
urite(1,1030) 'Contact ang. (deg) ', phi*180./pi
write(1,1000)

jsize - 50

jstep - 1

if (n/jaize .ne. 0) jstep - n/jsize

do 150, i-1, n, jstep
pctimpg - DSQRT(S(i))/(s2(i)*(1.-V(i)))
write(1,1050) i, x(i), 2(1), 0(1), Pctimpg, s2(i), V(i)
write(2,1090) 8(1), P(i)
write(4,1090) 8(1), pctimpg
urite(5,1090) 1(1), s2(i)
urite(6,1090) x(i), V(i)

if (i-jstep .ne. n) then
pctimpg ' DSQRT(S(n))/(82(n)*(1.-V(n)))
uncoudosm n. 3(a), 9(a). 0(a). Pctimpq: 32m). V(n)
write(2,1090) 3(a), P(n)
urite(3,1090) x(n), 0(a)
write(4,1090) 1(a), Pctimpg
write(5,1090) x(n), s2(n)
write(6,1090) x(n), V(n)
else
endif

186

write(1,1060) 'Pin ', Pin
write(l,1070) 'Pmax ', Pmax
write(l,1070) 'Xmax ', Xmax
write(1,1080) '0 (x - 0) ', 00
write(1,1070) '0 Lower Bound ', le
write(1,1070) '0 Upper Bound ', Qub
urite(1,1080) 'rorce ', forc
urite(1,1070) ’Load ', load
write(1,1080) 'Coat Thickness ', coat
urite(1,1080) 'tiber Vol tract ', Vt
write(1,1070) 'tiber Mass tract ', Mt
urite(1,1070) 'Resin Mass tract ', 1.0D+000 - Mf

close(unit-1)

close(unit-2)

close(unit-3)

close(unit-4)

close(unit-5)

close(unit-6)

999 format(/1x,'Tension Consolidation Iter.:',16)
1000 format(/3x,'Pressure Consolidation Iter.:',16)
1010 format(/T6,'iterQ',T10,'Q(0)',T31,'P(Q(0))'l)
1020 format(1x,18,2015.5)
1030 format(T15,324,T40,'-',T45,811.4)
1040 format(/T4,'Pt',T13,'x',T24,'P',T35,’Q’,T43,'§ IMPG',

8 756,'82',T68,'V'/)

1050 format(1x,IS,T10,P6.0,5P11.4)
1060 format(ll/T20,Al9,T40,'-',T45,F9.4)
1070 format(T20,Al9,T40,'-',T45,P9.4)
1080 format(/T20,h19,T40,'-',T45,r9.4)
1090 format(2t10.5)

end

***** sub: jacobian *****
subroutine AJRCOB(x,y,RJ)

implicit double precision (A-Z)

dimension RJ(2,2), y(2)

common c, DBLTAP, GAMMA, kappaP, kappaT, Patm
common Pin, 00, rf, $0, s20, theta, V0, va, Vt
external BETA, FV, Kt

V - FV(x)
call BETA(x,Bl,BZ,B3)

AJ(1,1) - o.o+ooo
AJ(1,2) - -c*B3/DSQRT(y(2))
AJ(2,1) - 2.*c*Kt(V)*(1.-V)
AJ(2,2) - o.n+ooo

return
end

181
***** sub: misc functs *****

subroutine BETA(x,BETA1,BETA2,BBTA3)

implicit double precision (A-Z)
external Ps2, WALL

s2 - Fs2(x)
H - WALL(x) + 1.D+000

BETAI - 1./3.*H**3 - 82**2*H + 2./3.*32**3
BETAZ - 1./2.*(ﬂ - :2)**2
BETA3 - -6./(H-82)**3

return
end

***** sub: modified jacobian *****
subroutine BJACOB(AJ,dx,BJ)

implicit double precision (A-Z)
dimension AJ(2,2), BJ(2,2)

a - 4.35866524D-001

BJ(1,1) - 1.D+000
BJ(1,2) - -1.*dx*a*AJ(1,2)
BJ(2,1) - -l.*dx*a*AJ(2,1)
BJ(2,2) - 1.D+000

return

and

***** sub: regula falsi - pressure consolidation *****
subroutine CONSLP(P,S,s2,V)

implicit double precision (A-Z)

common c, DBLTAP, GAMMA, kappaP, kappaT, Patm

common Pin, 00, rf, so, s20, theta, V0, Va, Vt

RBS(Vf) - (DSQRTWa/Vf) - 1.)": - kappaP/P*(DSQRT(Vf/V0) - 1.)
s1 - :2 - DSQRT(S)/(l. - V)

***** note: tol <- (DSQRT(Va/V) - 1.)**4 *****

tol - 1.D-010

aa - V0

bb - Va
cc - V0

10

10

188

if (DABS(RES(cc)) .ge. tol) then
cc - («mas (bb) -bb*RES (aa) ) / (RES (bb) -RES (aa))
if (RES(aa)*RBS(cc) .lt. 0.D+000) then
bb - cc
else
aa - cc
endif
goto 10
else
endif

V - cc
s2 - s1 + V0/V*(s20 - s1)

return
and

***** sub: regula falsi - tension consolidation *****
subroutine CONSLT(forc)

implicit double precision (A-Z)
common c, DBLTAP, GAMMA, kappaP, kappaT, Patm
common Pin, 00, rf, $0, s20, theta, V0, Vs, Vt

RES(Vf) ' (1./Vf - 1./Va)**2 - kappaT/forc*(1./Vt - 1./Vf)
***** note: tol <- (l./V - 1./Va)**2 *****
tol - 1.D-010

aa - Vt
bb - Va
cc - Vt

if (DABS(RBS(cc)) .ge. tol) then
cc - (aa*RBS(bb)-bb*RBS(aa))/(RBS(bb)-RBS(aa))
if (RBS(aa)*RBS(cc) .lt. 0.D+000) then
bb - cc
else
aa - cc
endif
goto 10
else
endif

V0 - cc
s20 - Vt/VO

return
end

10

189

***** sub: cubic spline fit of s2 and V *****
subroutine CSFIT(n,x,s2,V)

implicit integer (i,n)

implicit double precision (A-R,J-M,o-Z)

parameter (inc-50001)

dimension s2(inc), V(inc), x(inc)

common/CSPL/ xx(inc), 822(inc), VV(inc)
common/CSPLsZ/ Bs2(inc), Cs2(inc), Ds2(inc), ndats2
common/CSPLV/ BV(inc), CV(inc), DV(inc), ndatV
external SPLINE ,

do 10' 1-1, H
18(1) -x(i)
822(1) - s2(i)
VV(i) -V(i)

ndats2 - n
ndatV - n

call SPLINE(ndats2,xx,s22,Bs2,Cs2,Ds2)
call SPLINE(ndatV,xx,VV,BV,CV,DV)

return
end

tease sub: dfi/dx state
subroutine DFDX(dx,x,f,y,df)

implicit double precision (A-Z)

dimension df(2), f(2), y(2)

common c, DELTAP, GAMMA, kappaP, kappaT, Penn
common Pin, 00, rf, $0, s20, theta, V0, Va, Vt
external BETA, DWALL, Fs2, FV, Kt, MALL

V - BV(x)

s2 - Ps2(x)

H - WALL“) + 1.
call BETA(x,Bl,BZ,B3)

***** approx derivatives: ds2/dx, dV/dx, dB3/dx, th/dx *****

s22 - Ps2(x+dx)
V2 - FV(x+dx)

ds2 - (s22-s2)/dx

dV - (V2-V)/dx

th - (Kt(V2)-Kt(V))/dx

dB3 - BZ*B3**2*(DWALL(x) - ds2)

terml - f(1)/B3*dB3
term2 - -c*B3*(f(2)/DSQRT(y(2)) + DWALL(x) - ds2)

10

20

190

df(1) - terml + term2

terml - 2.*c*Kt(V)*(1.-V)*f(l)

term2 - -2.*c*Kt(V)*(y(1)-(Patm-GAMMA))*dV
term3 - 2.*c*(1.-V)*(y(l)-(PaterAMMA))*th
df(2) - terml + term2 + term3

return

end

***** fcn: derivative of wall boundary *****
sweat WEDGE tests

double precision function DWALL(x)

implicit double precision (A-Z)
common c, DELTAP, GAMMA, kappaP, kappaT, Patm
connon Pin, 00, rf, SO, s20, theta, V0, Va, Vt

DMALL - -c*DTAN(theta)

return
end

***** sub: max diff between full and half steps ***
subroutine ERROR(y1,y2,err,q)

implicit integer (i)

implicit double precision (A-R,J-2)
dimension err(2), to1(2), y1(2), y2(2)
column/ERROR] Etol

do 10, 1-1, 2
tol(i) ' y2(i)

do 20, i-l, 2
err(i) - y2(i) - y1(i)

q - DMAX1(DABS(err(1)/tol(1)1, DABS(err(2)/tol(2)))
q - q/Etol

return
end

***** sub: flowl (inlet flow bounds) *****
subroutine PLOW1(le,Qub)

implicit integer (i)

implicit double precision (A-H,J-Z)

parameter (inc-10000) '

dimension F1(inc+1), F2(inc+l), F3(inc+1)
common c, DELTAP, GAMMA, kappaP, kappaT, Patm
common Pin, 00, rf, SO, s20, theta, V0, Va, Vt
external BETA, Fs2, FV, SIMPSON

10

10

191

do 10, i-l, inc+1
x - (i-l)/DFLOAT(inc)
s2 - Ps2(x)
V - FV(x)

call BETA(X,81,82,B3)

F1(1) - 83

F2(1) 1./82

S (sZ*(1.-V))**2

P3(i) 33*(DSQRI(S) - DSQRT(SO))

call SIMPSON(inc,F1,P1int)
call SIMPSON(inc,F2,F2int)
call SIMPSON(inc,F3,F3int)

01b - l./F1int*(0.S*DELTAP/c - 3./2.*F2int)
Qub - 1./Plint*(0.5*DELTAP/c - 3./2.*F2int + F3int)

return
end

***** sub: integral property to find 00 *****
subroutine PLOW2(n,x,S,RES)

implicit integer (i,n)

implicit double precision (A-H,J-M,O-2)
parameter (inc-10000, isp1-50001)

dimension F1(inc+1), F2(inc+1), P3(inc+1)
dimension S(ispl), x(ispl)

dimension BS(isp1), CS(ispl), BS(ispl)

common c, DELTAP, GAMMA, kappaP, kappaT, Patm
common Pin, 00, rf, SO, s20, theta, V0, Va, Vt
external BETA, F82, PV, SPLINE, SEVAL, SIMPSON

call SPLINE(n,x,S,BS,CS,DS)

do 10, i-l, inc+1
x1 (i-1)/DFLOAT(inc)
SS - SEVAL(n,xl,x,S,BS,CS,DS)
s2 Fs2(x1)
V - l'V(x1)

call BETA(X1,BI,BZ,83)

F1(i) - 83

32(1) - 1./82

F3(i) C B3*(DSQRI(SS) - DSQRT(SO))

call SIMPSON(inc,F1,Flint)
call SIMPSON(inc,F2,F2int)
call SIMPSON(inc,P3,PBint)
RES - 00 - l./Plint*(0.S*DELTAP/c - 3./2.*F21nt + FBint)

return
end

10

20

192
***** sub: calc tow force and die load *****
subroutine FORCE(n,x,P,S,forc,load)

implicit integer (i,n)

implicit double precision (A-H,J-M,O-Z)

parameter (inc-10000, ispl-50001)

dimension dP(ispl), PP(ispl), F1(inc+1), F2(inc+1)
dimension P(isp1), S(isp1), x(ispl), P(2), y(2)
dimension BdP(ispl), CdP(ispl), DdP(ispl)
dimension BPP(ispl), CPP(isp1), DPP(ispl)

common c, DELTAP, GAMMA, kappaP, kappaT, Patm
common Pin, 00, rf, $0, s20, theta, V0, vs, Vt
external BETA, Ps2, SPLINE, SEVAL, SIMPSON, WALL

do 10, i-l, n
9(1) - P(i)
y(2) - S(i)
call FUNCT(x(i),y,f)
dP(i) - f(l)

call SPLINE(n,x,dP,BdP,CdP,DdP)
call SPLINE(n,x,P,BPP,CPP,DPP)

do 20, i-1, inc+1
x1 - (i-l)/DFLOAT(ino)
dex - SEVAL(n,x1,x,dP,BdP,CdP,DdP)
PP(i) - SEVAL(n,xl,x,P,BPP,CPP,DPP)
s2 - Ps2(x1)
ll - WALL(x1) + 1.

rl(i) - 1./(H - s2)
call BETA(xl,Bl,BZ,B3)
P2(i) - 1./c*dex*BZ*P1(1)

call SIMPSON(inc,Fl,P1int)
call SIMPSON(inc,F2,P2int)
call SIMPSON(inc,PP,load)

forc - Flint + FZint
load - load

return
and

****Q sub: eqn set *ittt
subroutine FUNCT(x,y,f)

implicit double precision (A-Z)

dimension f(2), y(2)

common c, DELTAP, GAMMA, kappaP, kappaT, PaUm
common Pin, 00, rf, SO, s20, theta, V0, va, Vt
external BETA, Ps2, FV, Kt, WALL

s2 - Ps2(x)

V - FV(x)

11 - MALL“) + 1.
call BETA(x,Bl,BZ,83)

193

Q - 00 + DSQRT(SO) - DSQRT(y(2))

{(1) - c*B3*(2.*Q - (8-3211
f(Z) - 2.*c*Kt(V)*(1.-V)‘(y(1) ' (Patm - GAMMA))

return
end

***** fcn: Kozeny constant *****
double precision function KOZ(V)

implicit double precision (A-Z)

Bl - 1.068136D-008
B2 - 2.5010250+001
B3 - 7.7435060+000

K02 - BI*DEXP(BZ*V) + B3

return
end

***** fcn: permeability *****

double precision function Kt(V)

implicit double precision (A-Z)

counon c, DELTAP, GAINA, kappaP, kappaT, Patm
common Pin, 00, rt, 80, s20, theta, V0, va, Vt
external xoz

Kt - (l. - V)**3/V**2*1./(4.*K02(V))*rf**2

return
end

***** sub: inverse of dim 2 matrix *****
subroutine MATINV(A,AINV)

implicit double precision (A-Z)
dimension A(2,2), AINV(2,2)

doth - A(1o1)*A(2o2) - A(1.2)*A(2,1)
detinv - 1./detA

Ainv(1,1) - detinv*A(2,2)
Ainv(l,2) - -detinv*A(1,2)
Ainv(2,1) - -detinv*A(2,1)
Ainv(2,2) - detinv*A(1,1)
return

end

10
20

10

20

194

***** sub: matrix mult *****
subroutine MATMULT(N,A,x,y)

implicit integer (i,j,N)
implicit double precision (A-n,K-M,o-2)
dimension A(N,N), x(N), y(N)

do 20, i-l, M
y(i) - 0.D+000
do 10, j-l, N

P(i) ' 3(1oj)*8(j) + V(i)

continue
continue

return
end

***** sub: 3rd order semi-implicit Runga-Kutta *****
subroutine RONGKUT(x,dx,y)

implicit integer (i)

implicit double precision (A-B,J-z)

dimension AJ(2,2), BJ(2,2), BJIMV(2,2)

dimension df(2), f(2), k1(2), k2(2), k3(2), y(2), Yl(2)
common c, DELTAP, GAMMA, kappaP, kappaT, Patm

common Pin, 00, rf, SO, s20, theta, V0, Va, Vt
external moon, Baacoa, omx, tuner, mum, MATMULT

a - 4.3586652400-001
b2 - 7.5000000000-001
b31 - -6.3020209OOD-001
b32 - -2.423378910D-001
Bl - 1.037609497D+000
32 - 8.349304840D-001

call AJACOB(x,y,AJ)
call BJACOB(AJ,dx,BJ)
call MATINV(BJ,BJINV)
call PUMCT (xoyo i)

call DPDX(dx,x,f,y,df)

do 10, i-l, 2
f(i) - f(i) + dx*a*df(i)

CALL HAMILT (2, BJINV, f, k1)

x1 - x + b2*dx
do 20, i-l, 2

y1(i) - y(i) + b2*k1(i)*dx

if (y1(i) .lt. 0.D+000) y1(i) - 0.D+000
continue

call PUNCT(x1,y1,f)

0000

30

40

50

10
20

195

do 30, i-l, 2
f(i) - f(i) + dx*a*df(i)

CALL MATMULT (2, BJINV, f, k2)

do 40, i-1, 2
terml - k1(i)*dx + a*df(i)*dx**2
term2 - k2(i)*dx + a*df(i)*dx**2
f(i) - b31*term1 + b32*term2

CALL HAMILT (2, BJINV, f, k3)

do 50, 1-1, 2
y(i) - y(i) + R1*k1(i)*dx + R2*k2(i)*dx + k3(i)
if (y(1) .lt. 9.0+000) y(1) - o.o+ooo

return
end

'**** fcn: determine $2 at x *****
double precision function Ps2(x)

implicit integer (i,j,n)

implicit double precision (A-H,K-M,O-z)

parameter (inc-50001)

common/CSPL/ xx(inc), s22(inc), VV(inc)
common/CSPLsZ/ Bs2(inc), Cs2(inc), Ds2(inc), ndats2
external SEVAL

Ps2 - SEVAL(ndats2,x,xx,s22,Bs2,Cs2,Ds2)

return
and

seval

from: G. Porsythe, M. Malcolm, and c. Moler
Computer Methods for Mathematical Computations
Prentice-Hall, 1977.

***** cubic spline interpolation - function *****
double precision function SEVAL(M,u,x,y,B,C,D)
implicit integer (i-k,N)
implicit double precision (A-H,L-M,o-Z)
dimension x(N), y(N), B(N), C(N), D(N)
DATA 1/1/
if (i .ge. n) i - 1

if (u .lt. x(i)) goto 10

if (u .le. x(i+1)) goto 30
i - 1

j - n+1
k - (1+J)/2

0000

0

3O

10

20

196

if (u .lt. x00) 1- k
if (u .ge. x(k)) 1 ' k
1: (j .gt. 1+1) goto 20
dx - u - 3(1)

SEVAL - y(i) + Dx*(B(i) + Dx*(C(i) + DX*D(i)))

return
end

***** fcn: non-linear spring *****
double precision function sigma(V)
implicit double precision (A-Z)

common c, DELTAP, GAMMA, kappaP, kappaT, Patm
common Pin, 00, rf, SO, s20, theta, V0, Va, Vt

sigma - kappaP*(DSQRT(V/V0) - lJ/(DSQRTWa/V) - 1.)“4

return
end

***** sub: integration - Simpson's rule *****
subroutine SIMPSON(inc,f,simp)

implicit integer (i)
implicit double precision (A-B,J-z)
dimension f(inc+1)

dx - 1./inc
simpO - f(l) + f(inc+1)
simpl - 0.D+000
simp2 - 0.D+000

60 10, 1-2, inc, 2
simpl - simpl + f(i)

d0 20' 1-3; inc-1' 2
simp2 - simp2 + f(i)

simp - dx/3.*(simp0 +4.*simp1 + 2.*simp2)

return
and

spline

from: G. Forsythe, M. Malcolm» and C. Moler
Computer Methods for Mathematical Computations
Prentice-Hall, 1977.

***** cubic spline interpolation *****

 

 

 

 

10

15

20

30

4O

50

197

subroutine SPLINE(N,x,y,B,C,D)

implicit integer (i,N)
implicit double precision (A-H,J-M,o-Z)
dimension x(N), y(N), B(N), C(N), D(N)

NMl - N-l
if (N .lt. 2) return
if (N .lt. 3) goto 50
0(1) - x(2) - 1(1)
C(2) - (V(Z)‘Y(1))/D(1)

do 10, i-2, NM1
D(i) - x(i+1) - x(i)

B(i) - 2.*(D(i-1) + D(i))
C(i+1) ' (Y(1+1) - V(i))/D(i)
C(i) - C(i+1) - C(i)

8(1) ' -D(1)

B(N) - -D(N-l)

C(1) - 0.

C(N) - 0.

if (N .eq. 3) goto 15
C(1) - C(3)/(x(4)-x(2)) - C(2)/(x(3)-8(1))
C(N) - C(N-1)/(x(N)-X(N-2)) - C(N-2)/(X(N-1)-X(N-3))
C(1) - C(1)*D(1)**2/(X(4)-x(1))
C(N) - -C(N)*D(N-1)**2/(x(N)-x(N-3))

do 20, i-2, N
T - D(i-1)/B(i-1)
3(1) ' 8(1) - T*D(1-1)
C(1) - C(1) - T*C(i-l)

C(N) - C(N)/B(N)
do 30, ib - 1, NMl
i - N - ib
C(1) - (C(1) - D(i)*C(i+1))/B(1)

B(N) - (y(N)-y(NH1))/D(NM1) + D(NM1)*(C(NM1) + 2.*C(N))'
do ‘0, i-l, NMl
3(1) ' (Y(1+1) ' Y(111/D(1) - D(i)*(C(i+1) + 2.*C(1))
D(i) ' (C(1+1) ' C(i))/D(i)
C(1) ' 3-*C(1)

C(N) - 3.*C(N)
D(N) - D(N-l)

return
3(1) - (P(Z) - Y(1))/(x(2) - x(1))
C(1) - 0.
0(1) - 0.
3(2) - 3(1)
6(2) - 0.
D(2) - 0.
return

end

198

***** sub: start integration (Euler method) *****
subroutine START(dx,x,P,S,Q)

implicit double precision (A-Z)

dimension f(2), y(2)

common c, DELTAP, GAMMA, kappaP, kappaT, Patm
common Pin, 00, rf, SO, s20, theta, V0, Va, Vt
external FUNCT

y(l) - P
y(2) - S

call FUNCT(x,y,f)

+ dx

+ dx*f(1)

+ dx*£(2)

+ DSQRT(SO) - DSQRT(S)

OM'UN
I I I I
OM'UN

***** sub: step size adjustment *****
subroutine STEP(q,dx)

implicit integer (i)
implicit double precision (A-H,J-z)

d1 - dx’DHIN1((4.*q)**(-.25),3.D+000)

return
end

***** fcn: determine V at x *****
double precision function EV(x)

implicit integer (i,j,n)

implicit double precision (A-H,R~M,o-Z)
parameter (inc-50001)

common/CSPL/ xx(inc), s22(inc), VV(inc)
comon/CSPLV/ BV(inc), CV(inc), BV(inc), ndatV
external SEVAL

EV - SEVAL(ndatV,x,xx,VV,BV,CV,DV)

return
end

199

c ***** fcn: wall boundary *****
C than WEDGE ﬁtter

double precision function WALL(x)

implicit double precision (A-Z)
common c, DELTAP, GAMMA, kappaP, kappaT, Patm
common Pin, 00, rf, SO, s20, theta, V0, Va, Vt

WALL - c*DTAN(theta)*(1. - x)

return
end

AA WE:

.OOOOOOD+000, 6.926950D+001

.899200D-002, 6.428700D-001, 9.070000D-001

.676200D-001, 1.0000000-003

.1451000-001, 3.000000D+001

.000000D+000, 0.000000D+000, 0.000000D+000, 0.000000D+000
.137000D+000, 1.8000000+000

.OOOOOOD-OOQ, 1.000000D-005, 1.000000D-010, 1.0000000-003

HHODGDN

theta
rf
kappaP
Ca

Pin
rhor
Ctol

c

Vt , Va

kappaT

phi

Pout , Patm , pct
rhof

Etol , Qtol, , Ttol

‘Q‘Q‘Q‘

LIST OF REFERENCES

LIST OF REFERENCES

Adams, K.L., B. Miller and L. Rebenfeld, "Forced In-Plane Flow of an Epoxy Resin
in Fibrous Networks,” Polym. Eng. Sci. 26, 1434-1441 (1986).

Adams, KL. and L. Rebenfeld, ”Permeability Characteristics of Multilayer Fiber
Reinforcements. Part 1: Experimental Observations,“ Polym. Eng. Sci. 12, 179-185
(l991a).

Adams, KL. and L. Rebenfeld, "Permeability Characteristics of Multilayer Fiber
Reinforcements. Part 11: Theoretical Model," Polym. Eng. Sci. 12, 186-190 (l991b).

Ahn, K.J., LC. Sefaris and LC. Berg, “Simultaneous Measurements of Permeability and
Capillary Pressure of Thermosetting Matrices in Woven Fabric Reinforcements, " Polym.
Campos. , 12, 146-152 (1991).

Ambari, A. , C. Deslouis and B. Tribollet, "Coil-Stretch Transition of Macromolecules
in Laminar Flow Around a Small Cylinder," 0mm. Eng. Commun. 29, 63-78 (1984).

Bascom, W.D. and LB. Romans, 'Microvoids in Glass-Resin Composites,” I & E C
Prod. Res. Dev. 7, 172-178 (1968).

Beavers, 6.8. and DD. Joseph, ”Boundary Conditions at a Naturally Permeable Wall,”
J. Fluid Mech. 30, 197-207 (1967).

Biller, P. , H.C. ﬁttinger and F. Petruccione, 'Consistantly Averaged Hydrodynamic
Interaction for Dumbbell Models in Elongational Flow, " J. Gram. Phys. 85, 1672-1675
(1986).

Bird, R.B. , W.E. Stewart and EN. Lightfoot, Transport Phenomena, Wiley, NY (1960).

Bird, R.B., C.F. Curtiss, O. Hassager and R.C. Armstrong, Dynamics of Polymeric
Liquids, Volume 2. Kinetic Iheory, Wiley, NY (1987).

Boger, D.V. , “A Highly Elastic Constant-Viscosity Fluid, " J. Non-Newt. Fluid Mech.
3, 87-91 (1977/1978).

200

201

Burdette, S.R., P.L Coates, R.C. Armstrong and R.A. Brown, ”Calculations of
Viscoelastic Flow through an Axisymmeuic Corrugated Tube Using the Etplicitly
Elliptic Momentum Equation Formulation (EEME)" J. Non-Newt. Fluid Mech. 33, 1-23
(1989).

Chai, M.S. and Y.I. Yeow, ”Modelling of Fluid M1 Using Multiple-Relaxation-Time
Constitutive Equations, " J. Non-Newt. Fluid Mech. 35 , 459-470 (1990).

Chilcott, MD. and LM. Rallison, ”Creeping Flow of Dilute Polymer Solutions Past
Cylinders and Spheres,” J. Non-Newt. Fluid Mech. 29, 381-432 (1988).

Chmielewski, C., K.L. Nichols and K. Jayaraman, “A Comparison of the Drag
Coefﬁcients of Spheres Translating in Com-Syrup-Based and PolybuteneBased Boger
Fluids," J. Non-Newt. Fluid Mech. 35, 37-49 (1990a).

Chmielewski, C., C.A. Petty and K. Jayaraman, 'Crossﬂow of Elastic Liquids through
Arrays of Cylinders,” J. Non-Newt. Fluid Mech. 35, 309-325 (1990b).

Chmielewski, C., K. Jayaraman and CA. Petty, ”A Theoretical and Experimental Study
of Resin Impregnation into Fiber Bundles Using a Wedge-Slit Die, " Presentation at ASC
3rd Tech. Conf., Seattle, Sept. 1988.

Cressely, R. and R. Hocquart, ”Birefringence D’écoulement Localisée Induite a L’arriere
D’obstacles," Optica Acta 27, 699-711 (1980).

Christiansen, R.L. and R.B. Bird, ”Dilute Solution Rheology: Experimental Results and
Finitely Extensible Non-Linear Plastic Dumbbell Theory, " J. Non-Newt. Fluid Mech. 3,
161-177 (1977/1978).

Christopher, R.H. and S. Middleman, ”Power Law Flow through a Packed Tube," I&EC
Fund. 4, 422-426 (1965).

Dabir, B. , ”Mean Velocity Measurements in a 3"-Hydrocyclone using Laser Doppler
Anemometry, " Ph.D. Dissertation, Dept. of Chem. Eng, Michigan State University
(1983).

Darcy, H.P.G., Les fontaines publiques de la villa de Dijon, Victor Dalmont, Paris
(1856).

Dean, M.M., Process Fluid Mechanics, Prentice-Hall, NJ (1980).

Dodson, A.G., P. Townsend and K. Walters, ”On the Flow of Newtonian and Non-
Newtonian Liquids through Corrugated Pipes” Rheol. Acta 10, 508-516 (1971).

202

Durst, F., R. Haas and B.U. Kacmaaar, ”Flows of Dilute Hydrolyzed Polyacrylamide
Solutionis in Porous Media under Various Solvent Conditions, " J. Appl. Polym. Sci. 26,
3125-3149 (1981).

Forsythe, G. M. Malcolm and C. Moler, Computer Methods ﬁrr Matheatical
Computations, Prentice-Hall, 1977 .

Gore, LM. and SP. Timoshenko, Mechanics of Materials, PWS Publishers, Boston
(1984).

Ghoniem, S.A. , “Extensional Flow of Polymer Solutions through Porous Media, " Rheol.
Acta 24, 588-595 (1985).

Gutowski, T.G., Z. Cai, S. Bauer, D. Boucher, J. Kingery and S. Wineman,
”Consolidation Experiments for Laminate Composites, " J. Campos. Mater. 21, 650-669

(1987).

Gutowski, T.G. , “A Resin Flow/Fiber Deformation Model for Composites, " SAMPE
Quarterly 16, 58—64 (1985).

Hudson, N.E., J. Ferguson and P. Maclde, 'The Measurement of the Elongational
Viscosity of Polymer Solutions Using a Viscometer of Novel Design, " Trans. Soc.
Rheol. 18, 541-562 (1974).

Inverarity, G. , ”Dynamic Wetting of Glass Fibre and Polymer Fibre, " British Polymer
J. 1, 245-251 (1969).

James, D.F. and DR. McLaren, “The Laminar Flow of Dilute Polymer Solutions
through Porous Media,” J. Fluid Mech. 70, 733-752 (1975).

James, D.F., N. Phan-Thien, M.M.K. Khan, A.N. Beris and S. Pilitsis, "Flow of Test
Fluid M1 in Corrugated Tubes" J. Non-Newt. Fluid Mech. 35, 405-412 (1990).

Kulicke, W.M. and R. Haas, ”Flow Behavior of Dilute Polyacrylamide Solutions through
Porous Media. 1. Influence of Chain Length, Concentration, and Thermodynamic Quality
of the Solvent," Ind. Eng. Chem. Fund. 23, 308-315 (1984).

Larson, R.B. and J .J .L. Higdon, ”Microscopic Flow Near the Surface of Two-
Dimensional Porous Media. Part 1. Axial Flow,” J. Fluid Mech. 166, 449-472 (1986).

Larson, R.B. and LLL. Higdon, ”Microscopic Flow Near the Surface of No-
Dimensional Porous Media. Part 2. Transverse Flow,” J. Fluid Mech. 178, 119-136

(1987).

203

Lawler, J .V., S.J. Muller, R.A. Brown and R.C. Armstrong, "Laser Dopler Velocimetry
Measurements of Velocity Fields and Transitions in Viscoelastic Fluids, " J. Nan-Newt.
Fluid Mech. 20, 51-92 (1986).

Lee, W.J., J.C. Sefaris and DC. Bonner, 'Prepreg Processing Science,” SAMPE
Quarterly 17, 58-68 (1986).

Lee, W., J. Seferis and L. Bravenec, 'Hot Melt Prepreg Processing of Advanced
Composites: A Comparison of Methods," SPE AMEC Tech. Papers 34, 1602-1607
(1988a).

Lee, W.J., LC. Sefaris and LC. Berg, ”Characterizing High Performance Composite
Processability with Dynamic Fiber Wettability Measurements, " Polym. Campos. 9, 36-41
(1988b).

Mackay, ME. and Ch.J.S. Petrie, ”Estimates of Apparent Elongational Viscosity Using
the Fibre Spinning and ”Pure“ Methods," Rheol. Acta. 28, 281-293 (1989).

Manero, O. and B. Mena, ”On the Slow Flow of Viscoelastic Liquids Past a Circular
Cylinder,” J. Nan-Newt. Fluid Mech. 9, 379-387 (1981).

Marshall, R.J. and A.B. Metzner, “Flow of Viscoelastic Fluids through Porous Media,”
Ind. Eng. arern. Fund. 6, 393-400 (1967).

McKinley, G.H., W.P. Raiford, R.A. Brown, and R.C. Armstrong, "Nonlinear
Dynamics of Viscoelastic Flow in Axisymmetric Abrupt Contractions, " J. Fluid Mech.
223, 411-456 (1991).

Mena, B. and B. Caswell, ”Slow Flow of an Elastic-Viscous Fluid Past an Immersed
Body,” Chem. Eng. J. 8, 125-134 (1974).

Michelsen, M.L., 'An Efﬁcient General Purpose Method for the Integration of Stiff
Ordinary Differential Equations,” A1015 J. 22, 594-597 (1976).

Middleman, S. , haldamentals of Polymer Processing, McGraw-I-Iill, NY (1977).

Muller, A.J., LA. Odell and LP. Tatham, ”Stagnation-Point Extensional Flow Behavior
of M1,“ J. Nan-Newt. Fluid Mech. 35, 231-250 (1990).

Nguyen, H. and D.V. Boger, ”The Kinematics and Stability of Die Entry Flows,“ J.
Nan-Newt. Fluid Mech. 5, 353-368 (1979).

Pilitsis, S. and A.N. Beris, "Calculations of Steady-State Viscoelastic low in an
Undulating Tube," J. Nan-Newt. Fluid Mech. 31, 231-287 (1989).

204

Pilitsis, S. and A.N. Beris, “Viscoelastic Flow in an Undulating Tube. Part II. Effects
of High Viscoelasticity, Large Amplitude of Undulation and Inertia, " J. Nan-Newt. Fluid
Mech. 39, 375-405 (1991).

Prilutski, G., R.K. Gupta, T. Sridhar and M.E. Ryan, ”Model Viscoelastic Liquids," J.
Nan-Newt. Fluid Mech. 12, 233-241 (1983).

Rallison, LM. and EL Hinch, "Do We Undestand the Physics in the Constitutive
Equation?, " J. Nan-Newt. Fluid Mech. 29, 37-55 (1988).

Richardson, S., “A Model for the Boundary Condition of a Porous Material. Part 2,” J.
Fluid Mech. 49, 327-336 (1971).

Sadowsld, T.J. and R.B. Bird, I'Non-Newtonian Flow through Porous Media. 1.
Theoretical," Trans. Sac. Rheol. 9, 243-250 (1965a).

Sadowsld, T.J. and R.B. Bird, ”Non-Newtonian Flow through Porous Media. 11.
Experimental," Trans. Sac. Rheol. 9, 251-271 (1965b).

Sangani, A.S. and A. Acrivos, ”Slow Flow Past Periodic Arrays of Cylinders with
Applieation to Heat Transfer,” Int. J. Multiphase Flow 3, 193-206 (1982).

Sangani, A.S. and C. Yao, ”Transport Processes in Random Arrays of Cylinders. II.
Viscous Flow,” Phys. Fluids 31, 2435-2444 (1988).

Sridhar, T., ”An Overview of the Project M1,“ J. Nan-New. Fluid Mech. 35, 85-92
(1990). .

Sridhar, T., R.K. Gupta, D.V. Boger and R. Binnington, "Steady Spinning of the
Oldroyd Fluid B. 11: Experimental Results,” J. Nan-Newt. Fluid Mech. 21, 115-126
(1986).

Taylor, G.I., “A Model for the Boundary Condition of a Porous Material. Part 1,” J.
Fluid Mech. 49, 319-326 (1971).

Villadsen, J. and M. Michelsen, Solution of Diﬂbrential Equation Models by Polynomial
Approximation, Prentice-Hall, NJ (1978).

Walters, K., 'The Seoond-Normal-Stress Difference Project,“ IUPAC MAORO 83,
Plenary and Invited Lectures Part 2, Bucharest, Romainia, 227-237 (1983).

Weaver, F., “Epoxy Adhesive Surface Energies Via the Pendant Drop Method,“ Tech.
Report AFWAL-TR-82-4l79, Air Force Wright Aeronautical Laboratories, OH, Dec.
1982.

205

Yeh, Y. and H.Z. Cummins, “localized Fluid Flow Measurements with an He-Ne Laser
Spectrometer," Appl. Phys. Lett. 4, 176-178 (1964).

“‘instillnu“