Skip to main content

Full text of "Modeling and experiment of bicomponent two-layer blown film coextrusion"

See other formats


MODELING AND EXPERIMENT OF 
BICOMPONENT TWO-LAYER BLOWN FILM COEXTRUSION 



By 



KUN-SANGjYOON 



A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL 
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT 
OF THE REQUIREMENTS FOR THE DEGREE OF 

DOCTOR OF PHILOSOPHY 

UNIVERSITY OF FLORIDA 



1993 



ACKNOWLEDGMENTS 



First of all, I would like to express my sincere gratitude to Dr. Chang-Won 
Park for his support. Without his academic advice and encouragement, this study 
could not have been accomplished. 

I would like to thank Drs. A. L. Fricke, L. E. Johns, R. Narayanan, and W. 
Shyy for their time on the advisory committee. Special thanks are extended to 
Sunkyong Industries of Korea and Department of Chemical Engineering at the 
University of Florida for their financial supports. 

I would like to extend my thanks to Do-Young and Woei-Shyong for their 
friendship and caring. 

Finally, I would like to share the accomplishment with my wife, Kyungha, 
and my daughters, Heewon and Jaewon and thank for their understanding through 
the study. 



n 



TABLE OF CONTENTS 



page 

ACKNOWLEDGEMENTS ii 

KEY TO SYMBOLS v 

ABSTRACT vii 

CHAPTERS 

1 INTRODUCTION 1 

2 STEADY STATE ANALYSIS 8 

Introduction 8 

Governing Equations and Boundary Conditions 10 

Nondimensionalization 20 

Numerical Computation 23 

Results and Discussion 25 

Summary and Conclusion 36 

3 STABILITY OF NEWTONIAN SINGLE-LAYER FILM 37 

Introduction 37 

Governing Equations 42 

Linear Stability Analysis 44 

Computational Procedure 49 

Results and Discussion 51 

Summary and Conclusion 62 

4 STABILITY OF BICOMPONENT TWO-LAYER FILM 64 

Introduction 64 

Governing Equations and Nondimensionalization 65 

Linearized Perturbation Equations and Boundary Conditions 68 

Computational Procedure 74 

Results and Discussion 78 

Summary and Conclusion 84 



5 EXPERIMENT 



87 



Introduction 87 

Experimental 88 

Results and Discussion 91 

Summary and Conclusion 100 

6 SUMMARY AND CONCLUSION 102 

APPENDICES 

A COMPUTER PROGRAMS FOR STEADY STATE SOLUTION 104 

B COMPUTER PROGRAMS FOR STABILITY ANALYSIS 110 

REFERENCES 136 

BIOGRAPHICAL SKETCH 139 



IV 



KEY TO SYMBOLS 



a = Metric tensor 

B = Dimensionless inflation pressure 

BUR = Blow up ratio 

e = Rate of strain tensor 



= Dimensionless rate of strain tensor 



= Take-up force 
= Flow rate ratio 



H 



= Melt thickness 



FL = Freeze line height 



M 

Q 

R 

Rf 



= Dimensionless freeze line height 
= Shear viscosity ratio 
= Volumetric flow rate 
== Radius of the bubble 
= Radius of the bubble at the freeze line 



Rl,Rh = Radii of curvature of bubble 



Tz 

t 

t 



= Dimensionless bubble radius 
= Dimensionless take-up force 
= Time 

= Dimensionless time 
= Thickness reduction, t^ = w(0)/w(^ ) 
= Velocity 

= Dimensionless velocity 



V 



w = Dimensionless melt thickness 

z = Axial length coordinate 

z = Dimensionless axial length coordinate 

AP = Pressure difference across the film 

A = Relaxation time of the UCM fluid 

A = Dimensionless relaxation time of the UCM fluid 
H = Shear viscosity 

cr = Total stress tensor 

a = Dimensionless total stress tensor 

T = Extra stress tensor 

T = Dimensionless extra stress tensor 

= Eigenvalue 

^ = Orthogonal curvilinear coordinates 



VI 



Abstract of Dissertation Presented to the Graduate School 
of the University of Florida in Partial Fulfillment of the 
Requirements for the Degree of Doctor of Philosophy 

MODELING AND EXPERIMENT OF 
BICOMPONENT TWO-LAYER BLOWT^ FILM COEXTRUSION 

By 

KUN-SANG YOON 
May, 1993 

Chairperson: Dr. Chang-Won Park 
Major Department; Chemical Engineering 

The flow mechanics of a bicomponent two-layer blown film coextrusion is 
studied theoretically. As a first step for the modeling of this complex process, a 
simple system is adopted in which the two layers are a Newtonian and an upper- 
convected Maxwell fluid (UCM), respectively. The two fluids are chosen to 
investigate the relative influence of viscous and viscoelastic forces on the overall flow 
of the process. For a given total flow rate, the radius and the melt thickness profiles 
of the blown film at steady state are determined numerically for various values of the 
flow rate ratio of the two fluids. When the relaxation time of the UCM layer is small, 
the flow mechanics is not much different from that of a Newtonian single-layer flow. 
With increasing relaxation time, the viscoelasticity effect of the UCM layer becomes 
increasingly pronounced and eventually dominates the bubble dynamics even though 
its layer thickness may be smaller than that of the Newtonian layer. 

The stability of the two-layer film blowing process is also investigated by the 
method of linear stability analysis. As a basis for the two-layer system, the stability of 



Vll 



a Newtonian single layer case is first considered for various processing conditions. 
Despite the simplicity of the model, the theoretical predictions appear to be in a 
qualitative agreement with many existing experimental observations. In the two-layer 
system the most significant effect of the viscoelastic Maxwell layer is that the 
maximum blow-up ratio for a stable operation is increased compared with that for the 
Newtonian single-layer case. The maximum thickness reduction for a stable 
operation, on the other hand, can be either smaller or larger than that for the 
Newtonian single-layer case depending on the magnitude of the viscoelasticity of the 
Maxwell fluid. Thus, if a material with moderately large viscoelasticity is chosen to 
be coextruded with a Newtonian fluid, the bubble stability can be maintained at a 
higher blow up ratio and a higher thickness reduction than possible with either 
materials alone. 

Bicomponent two-layer coextrusion experiments are also conducted to 
investigate the stability and to study the film physical properties influenced by the 
coextrusion. While the flow characteristics of the coextrusion process may be studied 
through a process modeling as presented in this thesis, the study on the film physical 
property is strictly empirical at the current stage of knowledge. A low-density 
polyethylene and a linear low-density polyethylene are used as the inner and the outer 
layer, respectively, for the experimental study. The results indicate that the addition 
of a viscoelastic low-density polyethylene layer tends to improve the bubble stability 
of the linear low-density polyethylene single-layer flow in accordance with the 
theoretical prediction. An adverse influence of low-density polyethylene, however, is 
observed in film physical properties in that the tensile strength of the linear low- 
density polyethylene is deteriorated by the presence of the low-density polyethylene 
layer. 



vni 



CHAPTER I 
INTRODUCTION 



The blown film extrusion is a process to fabricate thin plastic films with biaxial 
orientation by stretching the polymer melt in axial and radial directions 
simultaneously. Figure 1-1 shows a schematic of the blown film extrusion process. 
Polymer melt exiting an annular die is drawn upward to form a tube-like shape. Air is 
then introduced to the polymer tube through a hole in the middle of the die, which 
pressurizes the polymer tube thus inflating it to a bubble shape. In order to cool 
molten polymer and form a solidified film, air is supplied through an air ring which is 
located slightly above the die. The air ring distributes the air uniformly around the 
circumference of the bubble. Due to the cooling air, the bubble is solidified at a 
certain distance from the die, which is called a freeze line height. Between the die exit 
and the freeze line the polymer is in a molten state and all deformations (i.e., axial and 
radial stretching) occur in this region. Since the polymer is in a solid state beyond the 
freeze line, the bubble shape does not change any further although its temperature 
may still decrease due to cooling air. The solidified bubble then passes through nip 
rolls forming a flat double-layer film to be wound up for a post-conversion process. 
The nip rolls provide air seal so that the pressure difference between the inside and 
the outside of the bubble is held constant. 

The blown film extrusion is the most popular process for the manufacture of 
plastic films. In particular, a very large amount of polyethylene is consumed annually 
by the process to make thin films and sheets for numerous applications. Due to its 
importance in the plastics industry many investigations have been devoted to the 
mathematical modeling of the process as well as to studying the physical properties of 



1 



9 




Figure 1-1 Schematic of a blown film extrusion process 

R^, ; die radius Hq ; die gap 

Rf ; bubble radius H ; film thickness 
FL ; freeze line height 



3 



* 

the film as affected by various processing conditions. In a series of papers, 

Pearson and Petrie applied a fluid-mechanical analysis and developed for the first time 
a mathematical model for the process.**^ Their study was for a single-layer extrusion 
of an isothermal Newtonian fluid and numerous investigations followed to include the 
effects of heat transfer, gravity, inertia, and the viscoelasticity of the polymer melts. 

7 , 12-14 More recently, the influences of extrudate swell near the die exit as well as the 
transition from liquid-like to solid-like behavior near the freeze line have been also 
considered in modeling the blown film extrusion process. 

While the details of the flow mechanics may differ slightly depending on the 
complexity of the models that have been adopted, the overall features of the blown 
film process may be well described by the Newtonian isothermal model despite its 
simplicity. Due to the non-Newtonian characteristics of polymer melts and the 
difficulties in describing them using a simple constitutive relation, however, 
quantitative agreement between theory and experiments is not expected in general. 

Besides the steady state analyses, a few studies were also devoted to the stability 
analysis of the process, They were all based on the method of linear stability 

analysis, and only the axisymmetric disturbances were considered. 

During the past decade the fast advances in extrusion technology has also 
induced a very rapid growth in coextrusion processes. Coextrusion refers to a 
process in which two or more different polymers are extruded together to form a 
composite article. The coextrusion process offers a lucrative way to manufacture 
multilayer films with specific end-use properties such as barrier properties or heat 
sealability which are superior to those of single-layer films. Furthermore, compared 
to the multi-step lamination or coating processes, the coextrusion process is simple 
and economical and it now has become an indispensable process in plastics industry 
for those reasons. Despite its prevalence, however, there exist only a few studies on 
the coextrusion process. In the current study, a simple mathematical model for the 



4 



complex process is developed to investigate its steady state operation and its 
instability characteristics. 

The blown film coextrusion process is very complicated since the coextruded 
film usually consists of more than two layers made of different materials. All of the 
constituent materials are non-Newtonian with different rheological properties and the 
process is non-isothermal. A mathematical model which agrees quantitatively with 
experimental observations is rather specific to the choice of the materials for the 
experiment. Furthermore, its development may be an extremely difficult task due to 
the uncertainty in the constitutive equations. Reasonable qualitative predictions, 
however, can be obtained using a simple model as presented in this study. As a first 
step to tackle this complicated process, many simplifying assumptions are made. First 
it is assumed that the film consists of two layer made of a Newtonian and an upper- 
convected Maxwell fluid (UCM), respectively. Secondly, the flow is assumed to be 
axisymmetric, steady and isothermal. It is further assumed that only viscous and 
viscoelastic forces are important thus neglecting the effects of gravity, inertia, surface 
and interfacial tensions. 

The two fluids are specifically chosen to investigate the relative influence of 
the viscous and viscoelastic forces on the flow mechanics of the process considering 
the situation in which the two layers are a linear low-density polyethylene (LLDPE) 
and a branched low-density polyethylene (LDPE), respectively. It is well known that 
there exist fundamental differences in the flow behaviors of LDPE and LLDPE. A 
LDPE, which typically has a larger relaxation time than a LLDPE, shows an extension 
hardening behavior and tends to have better processability. Due to the rheological 
property differences, each material shows different bubble (i.e., tubular shaped blown 
film) dynamics when it is extruded as a single-layer film. If these distinctly different 
materials are coextruded together, it is anticipated that the bicomponent bubble 



5 



(whether it is two-layered or multi-layered) will behave differently compared to the 
single-layer ones of each material. 

It is obvious that a Newtonian and an UCM model cannot properly represent 
the LLDPE and LDPE, respectively. The two simple rheological models, however, 
may contain the essential features of the respective material in a qualitative sense. The 
former is a purely viscous fluid, whereas the latter provides a first approximation to 
elastic phenomena in fluids. The UCM model contains two material constants X and |i 

which represent the relaxation time and shear viscosity of the material. The 
elongational viscosity of the UCM fluid becomes unbounded when the extension rate 
is MIX. Thus, if X is small, the unbounded elongational viscosity is encountered only 

when the extension rate is very large. Consequently, the material behaves as if it is a 
Newtonian fluid for moderate values of extension rate. If X is large, however, the 
singularity in the elongational viscosity is encountered at a small extension rate, 
representing a strong extension hardening behavior. Although the present model is 
based on many simplifying assumptions, the current study can be easily extended in 
principle to a system in which more realistic but complicated rheological models are 
used and the effects of heat transfer, gravity, inertia, surface and interfacial tensions 
are considered. 

Much of the framework for the steady state analysis is similar to that of 
Pearson and Petrie for a single-layer blown film process. Following the theory of 
thin membrane (or "thin-film approximation"), the macroscopic force balances in axial 
and radial directions (z- and r-directions in Figure 2-1) are formulated. These force 
balances are solved numerically along with the constitutive equations to determine the 
radius and the melt thickness profiles as a function of material properties and 
processing conditions. When the relaxation time of the Maxwell fluid layer is small, 
the shape of the bubble is very similar to that of a Newtonian single-layer flow, 
indicating the negligible influence of viscoelasticity. With increasing relaxation time. 



6 



the viscoelasticity effect of the UCM layer becomes more and more important and 
eventually dominates the bubble dynamics. Thus, the two-layer bubble takes on the 
shape of a single-layer Maxwell fluid even though the flow rate of the UCM layer is 
smaller than that of the Newtonian layer. 

The stability of the system to an infinitesimal axisymmetric disturbance is 
investigated by the method of linear stability analysis. The operating parameters that 
determine the blown film extrusion are (1) the internal bubble pressure relative to the 
ambient pressure(B), (2) the amount of air trapped in the bubble (A), (3) the applied 
axial force (TJ, and (4) the take-up velocity(Vf). Either B or A in combination with 
either T^ or Vf can be fixed as the input parameters. Thus, in principle, four different 
cases can be considered for the blown film operation; 

Case 1 ; fixed B and Vf 
Case 2 ; fixed B and T^ 

Case 3 ; fixed A and Vf 
Case 4 ; fixed A and T^ 

Among these, the third one in which the amount of inflating air A and the take-up 
velocity v^ are fixed is the most relevant case to the industrial operations. Thus, the 
major emphasis of the present stability analysis is on the third case. The first and the 
fourth cases are also considered briefly to compare with existing studies. 

The bubble stability of a Newtonian single-layer extrusion is considered prior 
to the bicomponent two-layer extrusion since it serves as a basis for the more 
complicated two-layer case. The stability of the Newtonian single-layer case was 
investigated previously by two other investigators. Both studies, however, had a 
very limited scope in terms of the range of control parameters. Furthermore, critical 
mistakes were made in their mathematical formulation resulting in erroneous 
predictions. The present study not only covers a broad range of parameters that is of 



7 



practical interest, but also provides results that are in qualitative agreement with 
experimental observations. 

The stability characteristics of the bicomponent two-layer extrusion are quite 
different from those of the Newtonian single-layer extrusion in many respects. The 
most significant effects of the viscoelastic Maxwell fluid layer is that the maximum 
blow up ratio for a stable operation is increased whereas the maximum film thickness 
reduction is decreased. If the Maxwell layer is chosen to have a moderate level of 
viscoelasticity, however, bubble stability can be maintained for a larger operating 

range than obtainable with either the Newtonian fluid or the UCM fluid indicating a 
kind of synergism. 

The present study also includes two-layer coextrusion experiments, in which a 
linear low-density polyethylene (LLDPE) and a low-density polyethylene (LDPE) are 
used as the inner and the outer layer, respectively. These materials are chosen since 
they are known to have distinctly different rheological properties as well as film 
physical properties. The objectives of the experiment are to investigate processing 
characteristics (especially, bubble stability) for various values of the flow rate ratio of 
the two materials and to compare the physical properties of the film with different 
layer structures. The results indicate that the addition of viscoelastic LDPE layer 
tends to improve the bubble stability compared to the LLDPE single-layer extrusion. 
However, an adverse influence of LDPE layer on film physical properties is observed 
suggesting that caution should be exercised in selecting the materials for coextrusion 



processes. 



CHAPTER 2 

STEADY STATE ANALYSIS 



2.1. Introduction 

A mathematical model for the blown film extrusion process was first 
introduced by Pearson and Petrie. Using a thin film approximation in which 
velocity and stress fields are uniform across the film thickness, radial and axial force 
balances were derived for an isothermal Newtonian fluid system. Their analysis was 
then extended by several investigators to include the non-Newtonian and non- 
isothermal aspects of the actual process. Petrie attempted to investigate the influence 
of viscoelasticity and temperature variation using the Upper-Convected Maxwell 
(UCM) fluid model. ^ Due to the apparent mathematical difficulties involved with the 
numerical technique, however, his attempt was not successful. In section 2.4, details 
of the mathematical difficulties confronted by Petrie are discussed and a simple 
manipulation to resolve the problem is presented. Han and Park studied the power- 
law model for a non-isothermal case.^ More recently, Luo and Tanner investigated 
successfully the viscoelasticity effect using UCM and Luonov models. The 
influence of extrudate swell at the die exit as well as the phase change of polymer 
from liquid to solid near the freeze line was also investigated. 

In all these investigations, a shooting method along with the fourth-order 
Runge-Kutta method was employed for the numerical integration of the two-point 
boundary value problem. Recently, Cain and Denn adopted a new approach which is 
referred to as Newman's band matrix method, In that method, the ordinary 



8 



9 



difTerential equations for the system are reduced to a set of difference equations and 
the resulting algebraic equations are solved by matrix inversion. Both isothermal and 
non-isothermal cases were considered using three different rheological models; 
Newtonian, UCM and Marrucci model. The result for the Newtonian fluid is shown 
to be in good agreement with that of Pearson and Petrie. Their computation covered 
a broader range of variables which was not included in any other computations. 

The details of the flow mechanics may differ slightly depending on the 
complexity of the models that have been adopted; the overall features of the blown 
film process is well described by the Newtonian isothermal model despite its 
simplicity. Due to the non-Newtonian characteristics of the polymer melts and the 
difficulties in describing them using a simple constitutive model, however, quantitative 
agreement between experiments and theoretical investigation is not expected in 
general. 

During the past decade, significant progress has been achieved in extrusion 
technology and the use of coextrusion process has gained phenomenal popularity as a 
result. Despite the growing prevalence of coextrusion process, little attention has 
been given to mathematical modeling of the blown film coextrusion process. In this 
chapter, a rigorous analysis for the bicomponent two-layer blown film coextrusion 
process is presented for the first time using a simple mathematical model in which a 
Newtonian and an UCM fluid constitute the two layers. In the following section, 
governing equations and boundary conditions are derived for the two-layer case and 
the results of numerical computation are described for the steady state flow. The 
current study also uses the shooting method along with the fourth-order Runge-Kutta 
method as was done by Petrie.^ It is also presented that the mathematical difficulties 
encountered by Petrie can be resolved by a simple modification of the equations. The 
present results indicate that the effect of the viscoelasticity of the Maxwell fluid layer 



10 



becomes more significant with increasing relaxation time and eventually dominates the 
bubble dynamics. 

2.2. Governing Equations and Boundary Conditions 

The schematic of a bicomponent two-layer blown film coextrusion is shown in 
Figure 2-1. The inner and the outer layers are assumed to be a Newtonian fluid and 
an upper-convected Maxwell fluid (UCM), respectively. The region of interest for the 
present analysis is from the die exit to the freeze line, and the flow is assumed to be 
steady and isothermal throughout this region. For simplicity it is further assumed 
that only viscous and viscoelastic forces are important and the effects of gravity, 
inertia, and surface tension are neglected. 

Since the melt thickness is very small compared to the radius of the bubble, the 
thin-film approximation is applicable in which the velocity and the stress fields are 
uniform across the melt thickness. Thus the flow is locally a uniform biaxial 
extensional flow and only the diagonal components of the rate of strain tensor are 
non-zero. 

2.2.1. Kinematics 

The rate of deformation tensor for the steady state flow is derived in this 
section following the work of Pearson and Petrie. In the orthogonal cur\'ilinear 
coordinate system described in Figure 2-2, the local rate of deformation tensor is 



given as 



11 




Figure 2-1 Schematic of bicomponent two-layer blown coextrusion 

(H(z) ; film thickness, R(z) ; bubble radius) 



12 




Figure 2-2 Orthogonal curvilinear coordinate system which is used to 

describe the kinematics of the bubble 



13 




Here ^ 2 > ^3 represent the flow direction, the direction normal to the film, and 

the circumferential direction, respectively. The three components of the velocity 
vector of the fluid in the moving coordinate system are V,, Vj, and Vj. 

The velocity component normal to the bubble surface, Vj, is equivalent to the 
rate of thinning of the moving film element which is described as 

V2=dH/dt (2-2) 

Since H is a function of only ^j. 



aH 9H dt 9H , an 

-I- 2L = + Y 

at a^, at at 'a^, 



(2-3) 



By definition, V, =a|,/at. The equation (2-2) may be viewed as the material 
derivative for the fluid element convected with the velocity Vj in the direction as 
described in equation (2-3). The time derivative in equation (2-3) may be dropped for 
the steady state flow. Since the film thickness H is very small and the origin of the 
coordinate system is fixed at the inner surface of the film, the derivative of Vj with 
respect to may be given as 



14 




dH 



(2-4) 



The circumferential component of the velocity Vj is the expansion rate of the 
bubble in the circumferential direction. Therefore, 



■t r ~ dR dR 

V, = 2jt — = 27tV, 



dt 






(2-5) 



Then, the rate of deformation in the circumferential direction is given as 



V 3 ^ V, dR 
2 tcR ~ R d^, 



(2-6) 



Here R denotes the radius of the bubble (Figure 2-1). 

The local coordinate system is now transformed into the laboratory fixed 
coordinate system using the following transformation rule; 



dt = ^-dz = Vi + R'Mz 

' COS0 



(2-7) 



The rate of deformation tensor is then described as follows for the steady state: 



1 dV 




V dz 
0 



0 






0 

H dz 
0 



0 

0 

R dz^ 



(2-8) 



15 



In equation (2-8), V represents the velocity in the flow direction and the subscript 1 
for V is dropped since Vj and Vj are described using H and R. For the transient 
state, the diagonal components of the rate of deformation tensor are given as 



3R 

1-1- R'^ at ''' Vi-i-R'^ 

I 9H V 

H at 

1 aR . R' V 

R at ^ R Vl-I-R'- 



( 2 - 9 ) 



(2-10) 

( 2 - 11 ) 



Here the prime represents the differentiation with respect to z. The transient rate of 
deformation tensor is later used for stability analysis and a rigorous derivation of it 
may be found in Yeow.^^ 



2.2.2. The Continuity Equation 

Time dependent continuity equation is derived as follows by setting a material 
balance for a small fluid element: 



|-(rhVi-hr'') = |-(rhv) 
at dz 



For a steady state, this continuity equation is reduced to 



(2-12) 



Q = 2tcRHV 



(2-13) 



16 



Using this steady state continuity equation the velocity V in equation (2-8) is 
expressed in terms of H and R to give 






e = 



27tRHVl + R'' 






H dz R dz 
0 



0 



0 

1 dH 



H dz 
0 



\ 



0 

0 

R dz 



( 2 - 14 ) 



J 



It is the final form for the rate of deformation tensor which is used for the steady state 
computation. 



2.2.3. Force Balances 

Near the die exit extrudate swell typically occurs and the uniform flow 
assumption is not valid in that region. The extrudate swell, however, is restricted to a 
very small region in the immediate vicinity of the die exit and its influence may be 
negligible unless the thickness reduction (the ratio of initial to final melt thickness) is 
very small. Thus, the extrudate swell region is neglected in the present analysis. 

This situation is similar to the melt spinning or slot cast process for which the 
influence of the extrudate swell on the overall mechanics of the flow was proven to be 

negligible. 24, 25 , 27 , 28 

Under the prescribed conditions the shape of the bubble is determined by the 
balance between the viscous/viscoelastic forces in the film and the applied forces (i.e., 
axial tension and inflation pressure). The force balances in the axial and the radial 
directions can be given as follows: 



2jiRcos0(H“a;’, + H‘a|, ) + 7 i(r^ - R' )Ap = F, 



(2-15) 



17 



(hx, + h'<t;,) (hx, + hx,) 

Ap = + 



R 



R 



(2-16) 



H 



Here the superscripts o and i denote the outer and the inner layer, respectively. H is 
the melt thickness and R the radius of the bubble, a is the total stress tensor and the 
subscripts 11 and 33 indicate the normal stress components in the flow direction and 
in the circumferential direction, respectively. The first term in the axial force balance 
(2-15) represents the stress in the film, the second term is the net force in the axial 
direction resulting from the internal pressure of the bubble, and is the take-up force 
applied at the freeze line. 0 is the angle between the tangential line to the bubble and 
the vertical axis z (Figure 2-1). In equation (2-16) R^ is the radius of the bubble at the 
freeze line z=FL and Ap represents the pressure difference across the film. Rl and R„ 
are the principal radii of curvature of the bubble in the ^[^j-plane (or rz-plane) and in 
the respectively (Figure 2-2). In the laboratory fixed cylindrical 

coordinate system Rl and Rh are given as 




sec^0 

d'R/dz^ 



where sec0 = 



Ry = Rsec0 

^l+(dR/dz)’ 



(2-17) 



The two terms in each parentheses of equation (2-16) represent the normal 
force contributions of each layers in the flow and in the circumferential directions. In 
writing these force balances it is assumed that the velocity is continuous at the 
interface thus indicating the good adhesion between the two layers. 



18 



2.2.4. Constitutive Equations 

For the present biaxial extensional flow, only the normal stresses are involved and 
the stress-strain relation for the UCM outer layer is given as follows: 






tii+X 



dx 



\ 






dt 



^i--2e X- 
^11 



J 



= 2lie, 



(2-18) 



Here d/dt is a material derivative and the subscript ii represents the three normal 
stress components in ^ 2 > ^3 directions. Under the thin film approximation and 
using equation (2-7), equation (2-18) is reduced to 






r 



\ 



9t.. Vt' 

11 I *^11 



\ 



+ 



^ Vl + R 



n 



-2e,.X: 



11 ’ll 



J 



= 2|ie, 



(2-19) 



The total normal stress component in the film thickness direction is balanced with 
the external pressure (or ambient pressure) which can be set to zero without loss of 
generality. Therefore, Xjj is equal to the isotropic pressure p and consequently, the 
total normal stresses in the 1 and 3 direction are given as 




( 2 - 20 ) 



Using equation (2-20), equation (2-19) becomes 






3a 

11 



-h 



Va' 

n 



. ^ Vi 



-l-R 



/2 



2e-Gy 2X22(6;; e 22 ) 



2^(61.-622) 



(2-21) 



19 



Here the subscript ii is for 11 and 33 only. For steady state analysis, the time 
dependent terms are dropped and the velocity term are eliminated using the steady 
state continuity equation. 

If X is set to zero in equations (2-19) and (2-21), the equations for a 
Newtonian inner layer is obtained as 




X, = 2 tie, 

= 2 ti(eij —^ 22 ) 



( 2 - 22 ) 

(2-23) 



2.2.5. Boundary Conditions 

At the die exit the radius and the die gap are given as and Hg. It is assumed 
that the film is solidified at the freeze line and no further deformation occurs beyond 
this point. Thus the boundary conditions at each end are given as 

At z=0; R=Ro, H=Ho (2-24) 

At z=FL; dR/dz=0 (2-25) 

In addition to these conditions we also need to specify the initial stress condition for 
the viscoelastic outer layer as 

At z — 0 , — O’li.o *^33 ~ ^ 33,0 ^22 ~ ^ 22,0 (2*26) 

The stress state at z=0 accounts for the strain history of the fluid to that point. 
Unfortunately, the initial stress condition is not a measurable quantity and should be 
assigned rather arbitrarily. The influence of initial stress condition, however, 
diminishes rather quickly over a very short distance from the die exit. Consequently, 
the overall bubble dynamics is insensitive to the numerical value of the initial stress 



20 



condition. Thus the arbitrary choice of the initial stress condition may be 

justified. 



2.3. Nondimensionalization 



Proper scalings for the bubble radius R and the film thickness H are the radius 
and the gap of the die and Hq, respectively. The axial distance is also scaled by R^, 
and the isotropic pressure and all the stress components are scaled by viscous stress 
2Vop'/Ro- Under these scaling the force balances and the diagonal components of the 
rate of strain tensor are given in a dimensionless form as 



2frw{a„ +(l-f)(e„-e 22 )/f} = Vl + (rT(T + Br^) 
where T = T - B(BUR)^ 



(2-27) 






_-2Br 033 + (l - f)(e33 - )/ f 

T + Br^ r{a„+(l-f)(e„-e33)/f} 



(2-28) 



'^22 



<2 



.rwVl + (r') 



> = Me 



'22 ''22 



'22 



I — 2X22(6- ^22)} ^22) 

[rwVl + (r')’ 

where ii=l 1 and 33 



(2-29) 



(2-30) 



^11 =- 



1 



rw - J \ + (T') 



^w' r'^ 

— +— 

k W T J 



(2-31) 



1 






^22 ~ 



rw-^l-(-(r') 



k w; 



(2-32) 



21 




( 2 - 33 ) 



Here a and t are the total and extra stress tensors for the outer layer. In their 
expressions the superscript o is dropped since the stress-strain relationship for the 
Newtonian layer has been directly incorporated into the force balances. Prime 
denotes the differentiation with respect to the dimensionless axial direction z. r and w 
are the dimensionless bubble radius and the total film thickness, respectively. Once w 
is determined, the thicknesses of each layer are calculated by the conservation of 
mass: 



w°=fw. 



w'=(l-f)w 



( 2 - 34 ) 



The dimensionless parameters which appear in the above dimensionless 
equations are defined as follows: 



Dimensionless relaxation time: 



Viscosity ratio: 




Flow rate ratio: 



Dimensionless draw force: 



Bubble pressure: 



f = 



Q 



out 



total 



T. = 






B= 






22 



Blow up ratio; 


w 

c 

II 

O 


Freeze line height 


Ro 


Bubble radius 


R 

r = — 

Ro 


Fluid velocity 


< 

II 

o< < 


Film thickness 


K tif 

II 


Rate of strain tensor 


Vo 


Total stress tensor 


5= 

2VoP 


Extra stress tensor 


"fRo 

2V,tlo 


Axial length 


II 

IN 



Here X is the relaxation time of the Maxwell fluid. Vq and Q are the average axial 
velocity at the die exit and the volumetric flow rate, respectively. 

Finally the dimensionless boundary conditions are given as 



At z =0; 


r=l, w=l 




(2-35) 


At z =£ ; 


r'=0 (where t =FL/Ro) 




(2-36) 


At z = 0; 


^11 ~ ^11,0 ^ss ~ ^33.0 


"^22 ” ^22,0 


(2-37) 



23 



For the present analysis it is assumed that a,, o, Gjjg andx^jo are 0.1, 0.1 and 0.0, 

respectively. As pointed out previously, the overall bubble mechanics is insensitive to 
the initial stress condition which is rather arbitrarily assigned. 

2.4. Numerical Computation 

Equations (2-27) to (2-37) are a set of nonlinear ordinary differential equations 
for a two-point boundary value problem. Despite its complexity this set of equations 
can be solved by a shooting method once the material properties (X and M) and the 
operating conditions (f, T and B) are specified. First the value of r' at z=0 is 
assumed. Using this assumed value of r'(0) along with the initial conditions, the value 
of w' is obtained from equations (2-27), (2-31) and (2-32). Using this w', equations 
(2-28), (2-29) and (2-30) give the values of r", and ct'j, respectively. Thus all the 

derivatives are determined explicitly and they can be integrated from the die exit to 
the freeze line by a finite difference method. The initial value of r' is adjusted until r' 
at the freeze line is matched to be zero. The fourth-order Runge-Kutta method is 
used for the integration. 

If the flow rate ratio (f) is set to either zero or unity, the present system is 
reduced to the case of a single-layer Newtonian fluid or a single-layer UCM fluid, 
respectively. The results for either cases are known by several investigators who used 
at least two different numerical methods. Thus the accuracy of our computation can 
be compared with the existing results for the two single-layer cases and our results are 
shown to be in good agreement with those. 

The shooting method was proven to be effective for a Newtonian single-layer 
case.2.3 If it is applied to a single-layer UCM fluid, on the other hand, some 
numerical difficulties are encountered. If f=l, equation (2-27) is reduced to 



24 



2rw<T„ =7^ + (r')^(T + Br^) (2-38) 

Unlike the two-layer problem, this equation does not give an explicit expression of w' 
in terms of r' and o’,,. Instead w’ and ^ are coupled in the equations (2-29) and (2- 

30). Therefore at least two boundary conditions at z=0 (i.e., r'(0) and w'(0)) must be 
assumed in order to initiate the integration. Apparently this requirement causes rather 
severe numerical difficulties as encountered by Petrie^ and Luo and Tanner. To 
circumvent such difficulties, Cain and Denn employed the Newman's band matrix 

method. 

For the two-layer problem, however, the shooting method does not present 
such difficulties as described already at the beginning of this section. Furthermore, 
even for the case of a single-layer UCM fluid (i.e., f=l), it is found that a simple 
modification of equation (2-38) can in fact eliminate the difficulty of assuming two 
boundary conditions in using the shooting method. If we differentiate equation (2-38) 
with respect to z and use equation (2-28) and (2-30) to eliminate o'., an explicit 
expression for w' can be obtained in terms of r' as 



w = 



w 



((7,1 +4 T 22 + 2/a ) 



r'r"(T + Br=) 



2rwVl4-r'2 



w 




(2-39) 



Therefore if this equation is used instead of equation (2-38), only r'(0) needs to be 
assumed to initiate the integration and consequently the numerical difficulty 
encountered by the previous investigators is removed. The present results for the f=l 
case(UCM single layer) obtained by this shooting method is shown to be in good 



25 



agreement with those of Cain and Denn determined by the Newman's band matrix 
method.*"* 

It should be also pointed out that the specification of a,, at r=0 (i.e., equation 
(2-37)) is not necessary for the f=l case. As can be seen from equation (2-38), a,|(0) 

is determined once the initial value of r' is specified. For the two-layer problem, 
however, it is necessary to specify the initial value of ct,,. The numerical code written 

for the steady state solution is given in the Appendix section. 

2.5. Results and Discussion 

The variables that should be specified for the computation are the material 
properties A, M and the operating conditions f, Tz, B, and I (freeze line height). 
Once these six variables are specified, the computation results in the bubble shape 
r(z), and the thickness profile w(z). Thus the blow-up ratio (BUR) and the thickness 
reduction (t^) are also determined as a result. In an actual industrial process, BUR 
and tj are specified a priori. Thus the applied axial force T^ and the bubble pressure B 
are adjusted to obtain the desired values of BUR and U The large number of 
variables involved in the process not only generate numerous data points which may 
be presented in a number of different ways but also make the physical interpretation of 
the data a formidable task. Since our major interest is in investigating the relative 
influence of the viscoelastic force of the UCM layer to the viscous force of the 
Newtonian layer, the shear viscosity ratio M is fixed as unity while the values of f and 
X are varied. Considering the typical operating conditions in industrial practice, the 
values of axial force T^ and the inflation pressure B are determined so that the blow- 
up ratio (BUR) is bounded between 1 and 4 and the thickness reduction (t^) is smaller 
than 100. The thickness reduction t^ is defined as the ratio of the initial melt thickness 
to the final thickness. In addition, the value of (. is fixed at 5 for the results given in 



Bubble Radius 



26 




Axial Distance 



Figure 2-3 Radius profile of a single-layer bubble (t^=40; BUR=2.3 

The value in the parenthesis is TJ 



27 



this chapter. This value has been also chosen based on the conditions in industrial 
practice. 

The results obtained under the prescribed conditions are presented in two 
different types of plots; (1) the radius, the film thickness and the velocity profiles of 
the bubble and (2) the constant inflation pressure or the constant f contours in a BUR- 
t^ plane. In Figure 2-3 the shape of the single-layer bubble is given as a function of 
the axial coordinate z for various values of A for a fixed value of BUR and t^. Each 
curve, of course, has different values of B and T^ depending on the value of A . When 
A is smaller than about 0.2, the shape of the bubble is very close to that of a 
Newtonian bubble (A=0), indicating that the viscoelasticity effect is not significant. 
The slight downward shift of the A =0.1 curve compared to the Newtonian profile 
may be related with the initial stress condition which has to be specified unless A =0.0. 
When A is increased beyond 0.2, however, the bubble starts to take on a very 
different shape in that the initial rate of increase of the bubble radius is much larger 
than that of the Newtonian bubble. This change in bubble shape may be the 
manifestation of the viscoelasticity effect. As indicated in the figure, the required 
axial force (the values given in the parenthesis in the figure) also increases with A 
indicating the growing influence of the viscoelastic force. 

In Figures 2-4 and 2-5 the film thickness and the velocity profiles are given as 
a function of z. It seems that the larger value of viscoelasticity typically results in a 
fast reduction of film thickness at the early stage of elongational flow (i.e., for a small 
value of z) so that a large velocity gradient is prevented at the latter stage of the flow 
(i.e., for a large value of z). These predicted behaviors are common to other 
elongational flows such as fiber spinning or film casting. 22 similar behaviors 

have been in fact observed experimentally. 

The same trend is also observed for two-layer bubbles for which f=0.5 (Figure 
2-6). The bubble shape for A =0.3 is again very different from the others and the 



Film Thickness 



28 




Axial Distance 



Figure 2-4 Film thickness profile of a single-layer bubble 

(tj=40; BUR=2.3; Corresponding radius profiles 
are given in Figure 2-3) 



Velocity 



29 




Axial Distance 



Figure 2-5 Velocity profile of a single-layer bubble (t,=40; 

BUR=2.3; Corresponding radius profiles are 
given in Figure 2-3) 



Bubble Radius 



30 




Axial Distance 



Figure 2-6 Radius profile of a two-layer bubble (f=0.5; 

BUR=1.9; The value in the parenthesis is TJ 



31 



required axial force is also much larger than those for smaller A. Thus we may 
conclude that when A is 0.3, the viscoelasticity effect of the UCM layer is so 
dominant that the bubble takes on the intrinsic shape of an UCM fluid even though 
one half of the film consists of a Newtonian fluid. It may be also anticipated that, if A 
is larger than 0.3, the elasticity effect will become more pronounced and dominate the 
bubble dynamics even at a smaller value of f. 

For a given value of the inflation pressure B, there appears to be a minimum 
value of the axial draw force in order for a steady state solution to exist. The 
lower bound of the minimum T^ is approximately four times as big as B, which 
corresponds to the value for a straight cylindrical bubble (i.e., r=l) of a Newtonian 
fluid. When T^ is greater than the minimum value, two steady state solutions are 
typically obtained; one with a small BUR (consequently a small t^) and the other with 
a large BUR. This multiplicity of the steady state solution was first observed by Cain 
and Denn for the single-layer cases of both Newtonian and UCM fluids, In case 
of the Newtonian single-layer flow, the lower branch solution always has a BUR 
smaller than one which is of little practical interest. While two steady solutions are 
also obtained for the current two-layer problem, only the solutions with the BUR 
greater than one are considered in the present study. 

In Figure 2-7, the constant f contours are plotted on BUR-t^ plane for B=0.3 
and A =0.1. The f=0 contour represents the Newtonian single-layer flow. Along with 
this contour the axial force T^ increases monotonically from about 1.2 to 4.0 as t^ is 
increased from about 10 to 100, indicating a larger draw force for higher stretching. 
At BUR=1.0, Tz =1.2 corresponds to the value 4B which has been mentioned 
previously. Unlike the f=0.0 contour, the f=l contour for a UCM single layer flow is 
not continuous but splits into two branches. The lower branch with t^<40 is actually a 
closed curve if the data points with BUR<1 are added to this plot, and has much 
smaller value of T^ compared to the upper branch with t^>40. Along the upper branch 



Blow-Up Ratio 



32 




Figure 2-7 Constant flow rate ratio contours of steady state 

solution (dimensionless relaxation time A =0.1; 
dimensionless inflation pressure B=0.3; dimensionless 
freeze line height f =5) 



33 



the applied axial force T,, decreases with increasing t^. This trend (i.e., smaller axial 

force for higher stretching) is apparently non-physical. Thus this branch of solution 
may be unstable. 

When f-0.2 or 0.5, the shape of the contours is very similar to that of the 
Newtonian fluid for the given value of A . The T^ values for them are also within the 
same range as the Newtonian curve, indicating the insignificance of the viscoelastic 
force. The f=0.8 contour, however, resembles that for the UCM fluid (f=1.0), 
indicating the strong influence of the viscoelasticity effect. 

In Figure 2-8, a similar plot is given for a larger value of A (i.e., B=0.3, 
A =0.3). While the general trend is the same as Figure 2-7, a pronounced change is 
observed for the f=0.5 curve. Unlike the A =0.1 case, the f=0.5 contour now 
resembles that for the UCM fluid. It appears, however, that the contour is not yet 
separated but is in the transition from a single- to double-branched one. From these 
figures it may be concluded that, as A is increased, the elasticity effect of the UCM 
layer becomes stronger and stronger and eventually dominates the bubble dynamics 
even though its thickness may be smaller than that of the Newtonian layer. This trend 
is summarized in Figure 2-9 for a fixed value of B=0.3 and f=0.5. 

When f>0, at least one branch of each contours terminates at a point where 
BUR=1 and T^ is infinity. If the single-branched contours are extended below BUR=1 
(which is not shown in the figure), they all come down initially but curve back sharply 
upward to approach their respective asymptotic points. In case of the double- 
branched contours, the upper branch approaches directly to this point as shown in 
Figure 2-7 or Figure 2-8. The asymptotic point which corresponds to the solution for 
r=l (i.e., a straight cylindrical bubble) with an infinitely large applied axial tension can 
,in fact, be determined analytically. Since the radius of the bubble is constant as unity 
from the die exit to the freeze line and the bubble radius is much larger than in the film 
thickness, the flow situation is similar to the two-dimensional uniaxial elongational 



Blow-Up Ratio 



34 




Figure 2-8 Constant flow rate ratio contours of steady state 

solutions (dimensionless relaxation time A=0.3; 
dimensionless inflation pressure B=0.3; dimensionless 
freeze line height t =5) 



% 



35 




Thickness Reduction 



Figure 2-9 Constant flow rate ratio contours of steady state 

solutions (flow rate ratio f=0.5; dimensionless 
inflation pressure B=0.3; dimensionless freeze line 
height i =5) 



36 



flow for a flat film casting process. Thus, the asymptotic analysis for a bicomponent 
two-layer slot cast coextrusion process by Park is directly applicable to the present 
situation. The asymptotic analysis predicts that the asymptotic points for the f=0.8 
and 1.0 contours in Figure 2-7 to be located at t^=51.9 and 50.9, respectively. 
Similarly the asymptotic point for the f=0.8 and 1.0 contours in Figure 2-8 are 
predicted to be located at tr=17.1 and 16.1. 

2.6. Summary and Conclusion 

In this chapter the steady state flow of a bicomponent two-layer blown film 
coextrusion has been analyzed theoretically. Using a simple model in which a 
Newtonian fluid and an upper-convected Maxwell fluid constituted the two layers, the 
relative influence of the viscoelastic and viscous forces on the bubble dynamics has 
been investigated. Based on the theory of thin membrane, macroscopic force balances 
have been set up and solved numerically to determine the bubble radius and film 
thickness profiles as a function of material properties and processing conditions. 

When the dimensionless relaxation time (or Deborah number) is small, the 
bubble dynamics is not much different from that of a Newtonian single-layer flow 
since the influence of the viscoelasticity of the UCM layer is small. With increasing 
relaxation time, however, the viscoelasticity effect becomes so strong that it 
eventually dominates the bubble dynamics even though the flow rate of the UCM layer 
may be smaller than that of the Newtonian layer. Thus, the two-layer bubble basically 
takes on the shape of a single-layer Maxwell fluid under those conditions. The 
present model is based on many simplifying assumptions. Nevertheless, the general 
qualitative trends should remain the same even for the more complicated systems in 
which the flow is not isothermal and multi-layered. 



CHAPTER 3 

STABILITY OF NEWTONIAN SINGLE-LAYER FILM 



3.1 ■ Introduction 

The stability of the blown film extrusion process to an infinitesimal disturbance 
is investigated for the Newtonian single-layer case in this chapter. This investigation 
not only serves as a basis for more complicated bicomponent two-layer problem that 
is discussed in Chapter 4, but also provide useful information, in its own right, on the 
stability of the single-layer blown film extrusion. Due to its complexity, only two 
theoretical investigations have been attempted to date^^’l'*’^^, whereas several more 
studies exist on experimental investigation of the blown film stability.^^"^^ In general 
the stability of the biaxial extensional flow is limited by both the blow-up ratio (BUR) 
and the film thickness reduction (t^). In case of a linear low density polyethylene 
(LLDPE), the bubble is stable only up to BUR=2.0 whereas U can be reached as high 
as 100. High density polyethylene (HDPE) which has typically a higher molecular 
weight than LLDPE, on the other hand, sustains its bubble stability up to BUR^.O. 
Beyond these stable operating limits, bubbles show a few different instability modes 
before breakup; (1) an axisymmetric periodic fluctuation of the bubble radius and (2) 
a helical bubble behavior. 

Unlike the linear polymers, the stability of branched polymers seems to be 
more significantly limited by the film thickness reduction (t^). In case of low density 
polyethylene (LDPE), the maximum attainable t^ is about 30 whereas the maximum 
BUR can be much larger than that of LLDPE's. The instability mode is also 



37 



38 



apparently dilTerent in that the films usually breakup without distinct oscillatory 
behaviors. These differences are usually attributed to the higher viscoelasticity of 
LDPE resulting from its highly branched molecular structure. 

Among several experimental investigations, the work of Minoshima and 
White^^ seems to be the most comprehensive one. Using an annulus die with 1.6 cm 
in diameter, they investigated eight different polyethylene's; two LDPE's, four HDPE's 
and two LLDPE's. The blow-up ratio (BUR) and the film thickness reduction (t^) 
were varied to 5 and about 50, respectively, for two different freeze line heights (5 cm 
and 12 cm). These two freeze line heights correspond to the dimensionless freeze line 
height (() of 6.3 and 15.1. In addition to the two different modes of instability 
mentioned above (i.e., axisymmetric and helical instability), they observed another 
behavior in which the freeze line height fluctuated up and down periodically while the 
diameter of the bubble remained unchanged. They described it as a meta-stable 
behavior. While they pointed out the difficulty of physical interpretation of the 
experimental results due to the lack of theoretical understanding, their observations 
may be summarized as follows: 

(1) LDPE was the most stable whereas LLDPE was the least stable one among the 
three materials investigated. 

(2) The meta-stable behavior was observed only with one LDPE. 

(3) Increasing freeze line height decreased the range of stable operation. 

(4) The axisymmetric fluctuation of the bubble diameter appeared to be the most 
prevalent mode of instability which had been also observed by Han and his co- 
workers^^>3® and Kanai and White^f. 

(5) The helical instability which was observed at a high BUR might be associated with 
the cooling air blown along the outside of the bubble to control the freeze line 
height. 

In addition, Minoshima and White pointed out that the only existing theoretical 



39 



study at that time by Yeow showed little agreement with experiment. In Figure 3-1 
and 3-2, their experimental data for an LLDPE are plotted on a BUR-t^ plane for two 
different freeze line heights. The data for LLDPE are specifically chosen in these 
plots since the extensional rheology of this material is close to a Newtonian fluid. 
Although LLDPE shows the shear thinning behavior of most polymers, its extensional 
viscosity is known to be nearly constant due to its linear molecular structure. Figure 
3-1 does not have sufficient amount of data. Nevertheless, it seems apparent from 
Figure 3-2 that there exists an upper limit in both BUR and t^ Interestingly, the 
present study predicts such instability behavior at least qualitatively. 

In investigating the bubble stability, we may consider several different cases in 
which two of the four operating conditions B, T^, Vf and A may be held constant. 
Here B is the bubble inflation pressure (pressure difference between the outside and 
the inside of bubble) and is the axial force needed to draw the polymer melt (take- 
up force). Vfis the axial velocity of the polymer melt at the freeze line height (i.e., 
take-up velocity) and A represents the amount of inflating air trapped inside the 
bubble. Any two of these four parameters may be held fixed in theoretical analysis, 
whereas in practical operations the amount of air A and the take-up velocity Vf are 
always fixed. Among many possible cases, Yeow considered only one case in which 
B and Vf were fixed. Cain and Denn, on the other hand, considered four different 
cases: (1) fixed B and Vf, (2) fixed B and T^, (3) fixed A and Vj- and (4) fixed A and T^. 

Both investigations were linear stability analysis for an infinitesimal 
axisymmetric disturbance only, although their numerical techniques were different 
from one another. Yeow used a Newton-Raphson search technique, which located 
each eigenvalue individually for a given steady state value of B and T^. Cain and 
Denn, on the other hand, discretized the domain of integration and converted the 
eigenvalue problem into a set of algebraic equations, thereby obtaining the 
eigenvalues through a matrix manipulation. For the same problem (i.e., for fixed B 



40 




Thickness Reduction 



Figure 3-1 Stability characteristics of an LLDPE (Data taken 

from Minoshima and White^^ ; 

rio(Pa s)=l. 31x10'*; density=0.9181; dimensionless 

freeze line height f =6.3) 



Blow— Up Ratio 



41 




Unstable 



J L 



1 



i 



0 



10 20 30 

Thickness Reduction 



40 



Figure 3-2 Stability characteristics of an LLDPE(Data taken 

from Minoshima and White^^ • 

T|()(Pa s)=1.31xl0‘*; density=0.9181; dimensionless 
freeze line height f =15) 



42 



and Vf), their results were very different from one another. Yeow predicted that the 
film blowing process is always stable unless the thickness reduction is as large as 10'*. 
According to Cain and Denn's calculation, however, the bubble is mostly unstable 
except a small operating range where both BUR and t^ are small. Cain and Denn 
attributed this discrepancy to the difference in mathematical techniques. As will be 
described in the following sections, the present numerical results are shown to be 
closer to Cain and Denn's than Yeow's. 

The stability analysis of Cain and Denn covers only a single contour for B=0.2 
on a BUR-tf plane and does not provide the neutral stability curve which divides the 
BUR-tf plane into stable and unstable regions. In the following section, the stability 
analysis is repeated in order to provide a complete picture of the neutral stability 
curve within a typical operating range of BUR and t,. Furthermore, during this 
analysis a few critical mistakes have been found in formulating linearized equations in 
the work of Cain and Denn. 

The major emphasis of the present study is on the third case in which A and Vf 
are fixed. As pointed out previously this is the case relevant to the industrial 
operation. For the purpose of comparison the first and the fourth cases (i.e., (1) fixed 
B and Vfand (4) fixed A and T^,) are also considered. 

3.2. Governing Equations 

The time dependent continuity equation, the axial and radial force balances and 
the diagonal elements of the rate of strain tensor are given as follows; 

d d 

— (2;tRHV) = — — (2;rRHsec0) (3-1) 

oz at 



43 



— (2;rRH(T,,cos6- ;iR^AP) = 0 
dz 


(3-2) 


( rc rt \ 

AP = H 


(3-3) 


R' aR' V' 

^'■"l+R- at + 


(3-4) 


_ 1 an H' V 
Hat"^HVn-R'' 


(3-5) 


_ 1 aR R' V 

R at R Vl-hR'^ 


(3-6) 



The definitions of variables and the underlying assumptions in deriving these 
equations have been described in Chapter 2. In equation (3-2) differential form of the 
axial force balance is given since the applied axial force may fluctuate with time in 
a time-dependent transient flow. In case of the Newtonian single-layer flow the 
normal stresses can be expressed explicitly in terms of the rate of strain tensor. The 
dimensionless form of thus converted equations are given as 



^Jl+T 



/2 









<?w (9rA t ' dr' / / / « 

r-r:^-Hw— -t-rw-- 7 =====— ^-t-rvv'v +rwv-f-rwv = 0 

dt dtj Vl + r'^ ^ 



(3-7) 



2rr'w ar' 


2r 


aw 


(i/l + r'=f 


Vl + r'* 


at 



2rwv' 2rw'v „ , 

-Br^ 






= 0 (3-8) 



44. 




/ / / \ 
wv r w 




r V r w J 



(3-9) 



where the dimensionless time t is defined as t = tV(,/Ro, and prime denotes the 
differentiation with respect to z. As described in Chapter 2, r, w, and v denote the 
dimensionless bubble radius, film thickness and the velocity in the flow direction, 
respectively. 

In Cain and Denn's formulation the second term of the continuity equation 
(rwT '/Vl+(r')' )(*'/*■) was written erroneously as + Similar 

mistake was also made in the first normal component of the strain rate tensor, 
equation (3-4). Consequently, not only the continuity equation but also the force 
balances are erroneous. It is rather puzzling, however, that their results are often in 
qualitative agreement with the present results despite the obvious critical mistake. 



3.3.1. Linearized Perturbation Equations 

The standard procedure of linear stability analysis is applied to the current 
system, in which all the variables are expressed as a summation of steady state 
solution and a time-dependent small perturbation. 



3.3. Linear Stability Analysis 



r(zJ) = r„(z)-f-r'(z)e“' 
w( z, t) = w^ ( z) + w‘ ( z)e“' 
v(z,t) = v„(z)-Hv*(z)e“‘ 



(3-10) 

(3-11) 

(3-12) 



45 



Here the subscript ss indicates steady state, r*,w* and v* are the perturbed quantities 
and Q is a complex eigenvalue that accounts for the growth rate of the perturbation in 
the linear regime. Since only axisymmetric perturbation is considered, the variables 
do not depend on the azimuthal angle. In a transient state the internal bubble pressure 
B and the applied axial force can also fluctuate with time. Therefore, 



Equations (3-10) through (3-14) are substituted into the governing equations (3-7) 

through (3-9) and linearized about the small perturbation to give the following 
perturbation equations; 




(3-13) 

(3-14) 




(3-15) 



46 



2r„r' w 



ss ss ss 



(V'+(0')’ 



■n(r‘) 



2r 

ss 






£2w'+{2v,}n(w;) 






r(w V ) 4 - 2 Br 

,.2l' ss w ’ss ss/ ^*ss 



»■’ +{2(w„,fV;_f - }r; 



+ 



4r r' 



SS ss 



{‘+( 0 '} 



(w„vl, - wlv 



2^ ssss ”ss ss 



)(r-) 






2r V 



+ 






2r w' 



ss BS 



ss ss 






+{-2r,.twlr}v; +• 



.l + (r')=l 
-2r»w„ 



(w-) +{-2r^,fV^,f}(w;) 



•(''■) +{2ww^,f)(v;) 



(3-16) 




(3-17) 



47 



In equation (3-16) the subscript f indicates the variables evaluated at the freeze line 
height. This set of homogeneous linear differential equations may be solved by 
several different numerical methods to determine the eigenvalue once appropriate 
boundary conditions are specified. 

3.3.2. Boundary Conditions 

Part of the boundary conditions are obtained in a straightforward manner from 
equations (2-35) and (2-36) as 

r* = w* = v* = 0 at z = 0 (3-18) 

nr/ =-(rf*)'Vf at z = / (3-19) 

As mentioned previously, three different cases are considered; (1) fixed A and v^-, (2) 
fixed A and and (3) fixed B and Vf. Thus, appropriate conditions that describe 
fixed B, fixed T^ and fixed A needs to be given in addition to the above mentioned 
conditions. From equations (3-13) and (3-14), it is obvious that fixed B and fixed T 

z 

are equivalent to B*=0 and Tj,*=0, respectively. The description for fixed A, 
however, is not as obvious and requires a careful consideration. Since the absolute 
pressure inside the bubble is low, the amount of air trapped in the bubble can be 
expressed using the ideal gas law as 

(B + P)OT/L = nRT = A (3-20) 

Here (B+P) is the absolute pressure inside the bubble, P is the ambient external 
pressure, and L is the total height of the bubble from the die exit to the nip roll. Thus, 
lUfL represents the total volume of the bubble assuming that the height from the 

freeze line to the nip roll is much larger than the freeze line height. From the ideal gas 



48 



law n, R and T are the number of moles of the air, gas constant and the absolute 
temperature. The constraint that the amount of air A is constant is obtained by 
substituting equation (3-10) and (3-13) into equation (3-20) and linearizing it about 
the perturbation quantities; 



B-(qO + r;(2qP) = 0 



(3-21) 



Here it was assumed that B«P. This assumption is plausible since the internal bubble 
pressure is typically in the order of 10'^ atm in actual operations. 

Therefore, the appropriate boundary conditions for the three cases can be 
described as follows. 



Case I. 


fixed A and Vf ; 


B‘(r/) + r;(2rfP) = 0 


v;=o 


Case 2. 


fixed A and T^: 


B-(r,^) + r;(2qP) = 0 


o 

II 


Case 3. 


fixed B and v^; 


o 

II 

h 


< 

II 

o 



The conditions given in equations (3-18) and (3-19) are common for all these three 
cases. 

In the governing equations, T^ does not appear and it is necessary to express 
it in terms of other variables. This can be achieved as follows by perturbing the axial 
balance (equation (2-15)) at the freeze line height: 



(2r^)Qw* - 2(w^v; - w>^)r* +(2r^v;,)w“ +(-2r„v„)(w')* 

+ (-2r^<)v* -(2r^w„)(v')* -T; 



(3-22) 



49 



3.4. Computational Procedure 

The domain of our interest from the die exit to the freeze line height is 
discretized and the perturbation equations are reduced to a set of algebraic equations 
as follows; 



AiHx' =Bix' +Civ* 



A 2 QX' =B2X’ + C2V* 

where x’ =[r 2 ,r;,--",rj,w;,w;,----,w‘ 
and v' =[v;,v;,---,v;;.,,v>rT;]^ 




(3-23) 

(3-24) 



The subscripts for the vectors x’ and v* indicate the nodal points for which n 
corresponds to the freeze line height and 1 is for the origin of the axial coordinate. 

Ai and Bi are (2n-2)x(2n-2) matrices and Ci is a (2n-2)x(n-l) matrix whose 
components are determined from the finite differencing of equations (3-15), (3-16) 
and (3-19). A 2 and B 2 are (n-l)x(2n-2) matrices and C 2 is a (n-l)x(n-l) square 
matrix whose components are determined from equations (3-17) and (3-22). By 
eliminating v* between equation (3-23) and (3-24) the standard form for an algebraic 
problem is obtained as 



Ax*=Qx‘ (3-25) 

This equation is then solved using the subroutine GVLRG in IMSL package.^^ 

The first eigenvalue (i.e., the eigenvalue with the largest positive real part) 
depends on the spatial resolution or the number of discretization points (n) as 



50 



Table 3-1 Dependence of the eigenvalue with the largest real part on the number 

of discretization (The values for n— have been determined by the 
Richardson extrapolation method; FL=5, T^=2.74, B=0.45, 

BUR=1. 69, and t =22.29) 



n 


Or 


Oi 


41 


-0.3337 


11.4980 


51 


-0.3147 


11.5307 


61 


-0.3042 


11.5486 


71 


-0.2979 


11.5593 


81 


-0.2938 


11.5664 


n-^oo 


-0.2799 


11.5929 



51 



described in Table 3-1. While the limiting value of the eigenvalue for n->o« can be 
determined by the Richardson extrapolation method, most of the results presented in 
the following sections are obtained for n=51 to avoid excessive computation time. 
n=51 requires several inversions of a 100x100 matrix. 

3.5. Results and Discussion 

The results of the computation are plotted on blow-up ratio (BUR) vs. 
thickness reduction (t,.) plane. Figure 3-3 shows the stability of the case for fixed B 
and vf , which is the case considered by Yeow as well as Cain. The solid lines are the 
steady state constant B contours. The first and the second numbers in parentheses are 
the real part and the imaginary part of the eigenvalue with the largest real part, 
respectively. Contrary to the prediction of Yeow most of the typical operating region 

I 

is shown to be unstable to an infinitesimal axisymmetric disturbance except a small 
region where both BUR and C are small. In case of Cain's analysis only one contour 
for B=0.2 was investigated and the neutrally stable point was calculated to be at BUR 
=2.1 and t=4. This result appears to be quite close to the prediction of the present 
analysis although this agreement is somewhat puzzling since the Cain's formulation 
contains critical mistakes as pointed out previously. The imaginary part of all 
calculated eigenvalues are zero, indicating that the small perturbation either decays or 
grows without oscillation. 

The stability behavior of the bubble at fixed A and T^ is shown in Figure 3-4. 
The solid lines are the same constant B contours as in Figure 3-3, and the dashed line 
indicates the neutral stability curve above which the bubble is unstable. This figure 
indicates that under the given conditions the bubble is mostly stable if BUR is not 
larger than about 3.5. The real eigenvalues near the neutral stability curve may 
indicate a catastrophic failure without oscillation if the BUR is larger than a certain 



52 




Thickness Reduction 



Figure 3-3 The first eigenvalue of a Newtonian fluid single-layer 

case at fixed B and Vf (The first and the second number 
in the parenthesis represent the real and the imaginary 
part of the eigenvalue; Dashed line is a neutral stability 
curve above which the bubble is unstable; dimensionless 
freeze line height I =5) 



53 




Figure 3-4 The first eigenvalue of a Newtonian fluid single-layer 

case at fixed A and (The first and the second number 
in the parenthesis represent the real and the imaginary 
part of the eigenvalue; Dashed line is a neutral stability 
curve above which the bubble is unstable; dimensionless 
freeze line height i =5; dimensionless external pressure 
P=100) 



54 



critical value which is close to 3.5. 

The external pressure P which appears in equation (3-21) needs to be specified 
to simulate a situation for fixed A. For the present computations P has been set to 
100 which is based on the typical conditions of actual blown film operations. The 
numerical results based on P=100, 1000 and 3000 indicate that the location of the 
neutral stability curve is not seriously affected by the numerical value of P. Along the 
B— 0.2 contour, two transition points are found. As thickness reduction ratio is 
increased, the first transition from a stable state to an unstable state occurs near 
BUR=3.3 and tj-=24, and the second transition occurs at much higher thickness 
reduction ratio near tp=167 which is not displayed in Figure 3-4 since it is way outside 
the range of practical interest. 

The case of fixed A and vf is the most practical operating condition and the 
results for this case are shown in Figure 3-5. The constant B contours are again the 
same as previous figure and the dashed line also indicates the neutral stability. 
Compared to the previous case (fixed A and T 2 ) the stable region is much smaller. In 
addition to the unstable region where BUR is greater than about 3.5, the stability of 
the bubble is also limited by the thickness reduction the critical value of which appears 
to be increasing with BUR. Unlike Figure 3-4 the upper unstable region, however, 

has eigenvalues with non-zero imaginary part, indicating a different mode of 
instability. 

The neutral stability curves given in Figure 3-5 are very similar to those of 
Figure 3-2 at least qualitatively. Considering the fact that the case of fixed A and vf is 
relevant to the actual blown film operations, this agreement between the present 
results and the experimental observations by Minoshima and White is very interesting. 

It should be pointed out, however, that the experimentally observed upper unstable 
region is characterized by helical instability whereas the present analysis is for an 
axisymmetric instability. 



Blow-Up Ratio 



55 




Thickness Reduction 



Figure 3-5 The first eigenvalue of a Newtonian fluid single-layer 

case at fixed A and Vf (The first and the second number 
in the parenthesis represent the real and the imaginary 
part of the eigenvalue; Dashed line is a neutral stability 
curve; The bubble is stable in the region between the 
upper and the lower neutral stability curves; 
dimensionless freeze line height (. =5; dimensionless 
external pressure P=100) 



56 



Another interesting aspect of the present result is the critical value of t^ at 
BUR=1 which is shown to be 20.21. This value appears to be exactly the same as the 
critical draw ratio for the draw resonance in Newtonian fiber spinning or film 
casting. 2 ‘^. 27, 28 jp dimensionless bubble radius r is fixed at one, the governing 
equations for the axisymmetric blown film extrusion process become virtually the 
same as those for a two-dimensional slot cast process, and the same value of critical 
film thickness reduction should be obtained. In the present analysis, however, the 
bubble radius is allowed to fluctuate about its steady state value of one, yet the critical 
value of film thickness reduction is apparently the same as the critical draw ratio for a 
slot cast process. Furthermore, the imaginary part of the critical eigenvalue is also the 
same as that for a slot cast process, suggesting that the instabilities may be dictated by 
the same physical mechanism. 

In slot cast or fiber spinning process, an instability phenomenon called "draw 
resonance" is often encountered with linear polymer in which the film thickness or the 
fiber radius fluctuates periodically. While numerous experimental and theoretical 
studies have been conducted on this instability, the physical mechanism of this 
instability is not yet clearly understood. Hyun tried to explain this instability using 
kinematic wave theory.^'* He contends that there exists a throughput wave in the 
spin-lines which propagates from the die exit (or spinneret) to the take-up. The 
throughput wave resulting from disturbances travels with a unique wave velocity 
which is a function of draw ratio. If the residence time of the throughput wave is 
larger than one half of the spin-line residence time, the system is stable because the 
throughput wave (or disturbances) are washed out (or damped) due to the lack of 
time. Otherwise, the system is unstable. While Denn points out a number of errors in 
the mathematical manipulation of Hyun's analysis,^^ the basic idea behind Hyun's 
contention may be correct. The blown film process is much more complicated than 
the fiber spinning process due to additional parameters involved. Nevertheless, the 



57 



similarities in the governing equations for a special case of BUR=1 suggest that the 
instability mechanism may be the same. 

The influence of the freeze line height on the stability has been investigated in 
detail for the fixed A and v^ case since it is relevant to the industrial operation of 
blown film process. The dimensionless freeze line height C has been varied from 3 to 
12, and the results are shown in Figure 3-6 through Figure 3-9. In Figure 3-6, the 
neutral stability curves for £ =3 and £ =4 are plotted on a BUR-tr plane. The solid line 
is for £=3, and the dashed line is for £=A. As the freeze line height is increased from 
3 to 4, the area of stable region increases. When f=3, the bubble is found to be 
unstable when t^ is greater than about 30. When £ =4, however, the stable region is 
extended in the direction of higher BUR and higher t^ turning the Region I indicated 
in the figure into a stable region. The unstable region in Figure 3-6 is apparently 
divided into two regions for the given range of BUR and t^. Whether the upper and 
lower neutral stability curves will meet at a higher value of t^ is not clear since it is 
outside the range of the present analysis. 

The general trend of the neutral stability curve remains the same when the 
freeze line height is increased to 5 and 6 (Figure 3-7). Further increase in the freeze 
line height, however, results in a retraction of the stable region. When f=7, the 
neutral stability curve retracts back to where it was when £ =3, and further increase 
up to f =12 does not alter the curve significantly although the curves have a tendency 
to become more vertical as the freeze line height increases. The curve for f =15 is not 
given in Figure 3-8, since it almost coincides with the one for ^=12. It may be 
interesting to note that the upper neutral stability curve for £=5 is convex upward 
whereas that for £ =6 is convex downward. Not only the fact that the stable operating 
region is largest when £ =5, but also the change in the sign of the curvature may 
indicate the beginning of the retraction of the stable region with increasing £ . 

In Figure 3-9 three representative curves have been collected from Figure 3-6 



Blow-Up Ratio 



58 




Figure 3-6 Neutral stability curves for ^ =3 and 4 at fixed A and Vf 

(£ is dimensionless freeze line height; dimensionless 
external pressure P=100) 



Blow-Up Ratio 



59 




Thickness Reduction 



Figure 3-7 Neutral stability curves for £ =5 and 6 at fixed A and Vf 

{£ is dimensionless freeze line height; dimensionless 
external pressure P=100) 



60 




Thickness Reduction 



Figure 3-8 Neutral stability curves for E =7, 9 and 12 at fixed 

A and \^{E is dimensionless freeze line height; 
dimensionless external pressure P=100) 



61 




Thickness Reduction 



Figure 3-9 Neutral stability curves for £ =3, 5 and 7 at fixed A and Vf 

(£ is dimensionless freeze line height; dimensionless 
external pressure P=100) 



62 



through 3-8 to show the overall trend of bubble stability influenced by the freeze line 
height. For the investigated range of i (i.e., ^5) Region II is apparently always 

stable whereas Region I and III are always unstable. Region IV is shown to be stable 
only when 4<i <6. In industrial operations, it is desirable to have both BUR and t^ as 
large as possible to obtain a film with balanced physical properties. In that sense, the 
present result which suggests the existence of optimum freeze line height may have 
some practical implication. It may be also pointed out that the critical t^ is about 20 
when BUR=1 regardless of the freeze line height. As mentioned previously, this 
critical value may be closely related with that for the "draw resonance" in fiber 
spinning and slot cast processes. In these uniaxial extensional flow operations the 
critical draw ratio is also independent of the draw span which is equivalent to the 
freeze line height. It is also interesting to note that all the neutral stability curves 
nearly coincide when BUR<2. 

At the beginning of this chapter, the experimental study of Minoshima and 
White's was discussed, and their data for an LLDPE were plotted on a BUR-t^ plane 
to compare with the present theoretical results (Figure 3-1 and Figure 3-2). It is 
interesting to note the similarity between Figure 3-2 and Figure 3-9 in that the stable 
operating region is bounded by both BUR and t^. As pointed out previously, 
experimentally observed upper unstable region is characterized by helical instability 
whereas the present analysis is only for an axisymmetric instability. The angular 
dependence of perturbation may have to be included in the analysis to investigate the 
helical instability. 



3.6. Summary and Conclusion 



The stability of the blown film extrusion process for the Newtonian fluid single 
layer case has been investigated in this chapter. Among various cases that have been 



63 



considered, the case of fixed A (the amount of air trapped in the bubble) and Vj- (take- 
up velocity) is most relevant to the typical operations in industrial practice. Within 
the practical range of BUR and t^ (i.e., 1^UR<4 and t^<60), the present analysis 

predicts that the bubble stability is limited by both the blow up ratio (BUR) and the 
film thickness reduction (t^). This trend is apparently in qualitative agreement with 
experimental observations. Neither the mechanism of bubble instability nor its 
complex trend can be explained physically at this point. Nevertheless, the present 
study which predicted a realistic trend of bubble stability for the first time should 
provide useful information for future studies to refine the present simple model. 



CHAPTER 4 

STABILITY OF BICOMPONENT TWO-LAYER FILM 



4.1. Introduction 

The stability of a bicomponent two-layer blown film coextrusion is investigated 
in this chapter. In the previous chapter, the Newtonian single-layer case was 
investigated for various operating conditions. It was pointed out that among the 
various cases, the case of fixed A (amount of air trapped in the bubble) and Vf (take- 
up velocity) is most relevant to industrial practice, for which the stable operation is 
shown to be limited by both the blow-up ratio (BUR) and the film thickness reduction 
(tj). It was also pointed out that such stability behavior is in fact observed at least 
qualitatively with linear polymers such as LLDPE. In case of a LDPE which has a 
larger viscoelasticity than a LLDPE, on the other hand, experimental observations 
indicate that a higher BUR than obtainable with a LLDPE can be achieved whereas 
the maximum thickness reduction is generally smaller than that for a LLDPE. 

When these two materials are coextruded, the competing influence of the two 
fluids may significantly alter the stability behavior of the two-layer film. The objective 
of this chapter is to investigate the stability of a bicomponent two-layer film as 
affected by the rheological properties of two different materials. As discussed in 
Chapter 2, the present model consists of a Newtonian and an upper-convected 
Maxwell fluid as the inner and the outer layer, respectively. These purely viscous and 
viscoelastic constitutive models are chosen with LLDPE and LDPE in mind although 
neither models can exactly describe these materials. 



64 



65 



The mathematical approach for the two-layer film is in principle the same as 
that for the single-layer case of Chapter 3 although the added complexity is 
formidable. Only the case of fixed A and Vf are considered for the present analysis 
while the time constant for the UCM layer A and its flow rate fraction f are varied. 

The results indicate that the presence of the viscoelastic skin-layer removes the 
upper bound on BUR for a stable operation whereas the maximum film thickness 
reduction can be either smaller or larger than the value for a Newtonian single-layer 
case depending on the magnitude of A. The most interesting conclusion may be that 
at a low value of f (say f=0.2) the region of stable operation can be larger than those 
for a single-layer flow of either constituent materials. Thus, a kind of synergism is 
predicted for coextrusion in terms of bubble stability. 

4.2. Governing Equations and Nondimensionalization 

The time dependent continuity equation, the axial and the radial force balances 
for the two-layer flow are as follows: 



^(2nRHV) = -^(27tRHsec0) 
oz dt 

4(2;rRcose(HXi + - ;rR*Ap) = 0 

OZ 









(4-1) 

(4-2) 



(4-3) 



The definitions of variables have been given in Chapter 2. The continuity equation (4- 
1) is exactly the same as that for a single-layer flow, equation (3-1). In this equation, 
however, H is the total film thickness that is the summation of the inner and the outer 



66 



layer film thickness. The axial and the radial force balances (4-2) and (4-3) are also 
very similar to those of the single-layer flow, equations (3-2) and (3-3). The only 
difference is that the normal forces in the two-layer film are expressed as the 
summation of the inner and the outer layer contributions that are denoted by the 
subscripts i and o, respectively. Time derivatives do not appear explicitly in the force 
balances, but these equations are time dependent since the normal stress components 
vary with time; 









T^2 ^ 



dvl. 



+ 



V 



V 



dt Vl + R'^ 



2e T® 
^^22 ‘■22 



= 2|i®e 



22 



J 



(4-4) 



<7®i +A 



dal Vo® ^ _o o r \ 

+ , 2e,jO®i ^'^22i^ii ^ 22 ) 



^ Vl + R'^ 



= (e^, - e,,) 



'^22 '^^^22 

o', 



where ii=l 1 and 33 



(4-5) 

(4-6) 

(4-7) 



Equations (4-4) through (4-7) are the constitutive relations for the UCM and 
the Newtonian fluid. The diagonal components of the strain rate tensor are given in 
equations (3-4), (3-5) and (3-6). 

The final form of the governing equations for the two-layer flow is obtained by 
substituting the expressions for the strain rate tensor (equations (3-4), (3-5) and (3- 
6)) into the above equations. The dimensionless equations thus obtained are 



67 



VI 



/ 



4- 



V 



dw 

w — 

dt at 



r' 3 r' , , , 

+ rw --7======== + rw v-f r wv = 0 



y 



Vl + r'^ 



(4-8) 



Vi 



2frw [_ 

i^r" 



4- 



l~f 



/ 



r' ar' 1 aw 



4 - 



1 



/ 



V 



i4-r'^ at w at 






V 

V W 






I 



-Br 



= 0 (4-9) 



" ^ (VTi^) ^ 



r' 




[Vl + r'= 


r 

/ 



9w 



= B(l + r-) + f^2lL + (l-f)'^’'"''' 



Vi 



+ r 



/2 



1 + r 



/2 



-( 1-0 



/ 

wr V 



l + r 



f2 



(4-10) 



fZ^7l + r'^ _(l_f) 



r wv w V 

-^+(1-0— 



^(2Xr^„+2Xr't,.r-)|+l(-2XT„-l)^+(-X)fl 

l + r dt w dt dt 

__ ^va,', 2Xv'a,'i 2Xv'x,2 

Vl + r'^ Vl + r'^ Vl + r'^ 

2Xw'vt„ v' w'v 

-I — -j 

wVl + r'^ Vl + r'^ wVl + r'^ 



(4-11) 



-(2X033 +2XT22 + 1)^ - — (2XT22 + l)% + (-^)%^ 
r dt w dt dt 

- ?T -4- ^^^33 _ 2XrVa33 _ 2 Xr vTjj 

Vl + r rVl + r'^ ~rVl + r'^ 

_ 2Xw'vT22 ^ , w V 

wVl + r'^ rVl + r'^ wVl + r'^ 



(4-12) 



68 



w dt dt 

Xvt', 2XwVf,, w'v 

r: X H 

Vl + r'^ w-v/l + r'* wVl + r'^ 

The first three equations are the continuity equation, the axial and the radial force 
balances, respectively and the following three equations are the constitutive relations 
for the viscoelastic outer-layer (i.e., UCM fluid). Superscript o is dropped because all 
the stresses (cr,,, Tj,) in the dimensionless equations represent those of the 

outer layer. 



4.3. Linearized Perturbation Equations and Boundary Conditions 



The linearized disturbance equations are obtained by following the standard 
procedure of linear stability analysis. As in section 3.3, all variables are expressed as 
a summation of steady state solution and a time-dependent small perturbation 
(equation (3-10) through (3-14)). These variables are then substituted into the 

governing equations which are linearized about the perturbation to give the following 
equations; 





(wX+<v 
(r>„ + r„w;, 



«)r* +(v„w^)(r*)' + (rX + +( 

)v*+(r,,w^)(v*)' 







(4-14) 



69 



2(l-f) 



r„w„r„ 



(Vl + r'^ ) 



/ \ ♦ 



Q(r') 



+ 






» 



yTl 









Qw' +[2(l-f)r„f]Qw; 



-2f -2(l-f) - — . r “^ + 2Br 



Vi 



+ r 



/2 



ITTT^ 



SS 






.2fw^.ff^u,ss,f +2(l-f)(w„,fV;f - w;_fV„f)-2Br„,f]r; 



+ 



2 f [ 4(^1 



(yl^ + ^L^) 



3 



(V'+''«0 



r( 



w v' - w' V ) 

SS SS SS SS / 



/\ * 



(r') 



+ 



r„ O’, 1 „ r V ^ 

■7f. “ "■^ . -2(i-f) ^ 






f2 

SS 






W’ + [2fr^,f CT„^ f + 2(1 - f )r^f v^f ]w; 



+ 



r V 

2(l-f) “ “ 



(^) 



(w')*+[-2(l-f)r„fV^f](w;)’ 



+ 



r w 

2(l-f) " “ 






v*+[-2(l-f)r^fW;f]v; 



+ 



-2(l-f) 



r w 

SS SS 






(v')‘+[2(l-f)r„fW^f](v;)* 



r w 

_2f ss'Vs* 



■\/i+^_ 

+[r» -r«,f]B' 



(^n+[2fr«.fW„f](r;,f 



(4-15) 



70 



( 1-0 






.2 



ss 



nr' + 



-( 1-0 



/ // 



1 + C 




Q(r*)' 



+ 






( 1-0 



// 



V'+c 



\ 



V 



■\/i+C 



ss 



7 



Dw* = 






f ^ , ^ 

2r«w^ 

— — — — — Yv 

ss 



ss 



ss 






/ 



+ 



2BC 



f (w..v^ 




1 + r 



/2 



ss 



(i + Cj 



,2^ ' “ “ 



W V 

ss ss 



)-f 



'■ssW„(T33.ss 



-( 1-0 



SS 



ss ss 



ss 



• \/ 



(r) 



+ 



Vl + r/ 



./2 

ss 



1 + C 



/ / \ 

,W V — W V ) 
2 v'^ss'^ss '^ss'^ss/ 



• \ // 



(r ) 



+ 



-^2%+(i-f)-^-f^,/ri^-(i-f)% 

-y 1 + ^ ^s$ ^SS ^ss 



w 



+ 



-(1-f) 



CVss 

(l+'I') 



+(l-f) 



ss 



ss 



♦ \ / 



(w ) 



+ 



-(1 - f _ (I - f ) _L(r;^„_, - ) 



1 -f r 



ss 



ss 



V + 



( 1-0 



// 

r w 

ss ss 



1 + r 



/2 



ss 



• \/ 



(V^) 



+ 



+ 



// 

r IV, 

T ss ss 

^/l^ 



l + r'^]B' 







c\, + 






T 




L ss J 



33 



(4-16) 



71 



2 >.-!»-^( 5 „,+V„ + i) 



l + r 



n 



ss 



a(r*y 



+ 



1 /.T-_ 



w 



+ 1) 



ss 



+^-X]C!a,*, = 



ss 



(^/^^) 



■^^ss^ll.ss ^^^ss(^ll,ss "^22,ss) 



W>s,^22.ss 



+ 



W 



/ / M 

V' ^ssVss 
^ ss 



ss 



V 



w 






* \f 



(r) 



+ 



W V 

ss ss 



w 



ss 






/2 

ss 



+ 1 ) 



W 4- 



ss 



w 



ss 






/2 

ss 






(w-)' 



+ 



X 






+ 



w 



ss 



/2 

ss 



w 



ss 



VI77 



/2 

ss 



+ 1 ) 



V + 



-1 



.Vi 



+ r 



(2Xc,,^„ + 2 Xt22_^ + i) (v*)' 



ss 



+ 



W 

— V 

1-2X 



•A+n 



/2 

ss 





/ 


a,*,+ 


— V 

^ . .. ..“ - 


r^c\ 



(^u )' 



+ 



2X 



Vii^ 



/2 

ss 



w'v " 

-_y' 5 L_H 

ss 



V 



w 



ss 



'll 



(4-17) 



72 



1 i.^. 



( 2 X <733 ^ + 2 X,X 22 ,ss 



S5 



fir' 



+ 



1 /_T-_ 



w. 






ss 



nw' + [-X]Oa„„ = 



--4^r(2)i5,,,. + 2Xv„+l) 

/ssV^'^^ss 



-f 



\ ^ss^ss^33,ss 

’ (7T^)“ 



ss 






(2^CT33.«+2Xt22_„ + l) 






/2 A 



1 - 



ss 



V 



14-r 



/2 



ss / 



/ / 

r w V 

ss ss ss 






(23.t;j^ + l) 



w 



ss 



♦ \ / 



(r ) 



+ 



W V 

ss ss 



U' 



w:. a/ I + r 



n 



+ 1) 



SS 



w + 



ss 



Wss-\A + r«' 



(2^'C22.ss + 0 



(w'r 






(2^^33,55 2 XXjj.ss l) 



+ A 

+ C ^ssa/I+C 



w 



ss 



W 



ssaA^ss 



/2 



(2^'^22,ss + 0 



+ 



— r V 

1 _ 2 A- “ ” 



ss 



VTk 



^2 

ss 



<^ss + 



X 






/2 

ss 



iKY 



+ 



— V 

2X “ 



w„ r 






ss ss 



w 

\ '^ss 



‘ssy 



'22 



(4-18) 



73 



1 






( 2 -^ 'T’224. + 0 



Qw* +^-AjQT22 ~ 



-A 






SS S3 22,53 ^ SS SS SS 

U ^+ c ) "'.( v >+<0 






• \/ 



(r ) 



+ 



A 



x' 

‘'22.SS 



W 



(2-^T22,„ + 1 ) 



SS 



Vf 



,/2 



+ r„ w„ a/ 1 + r 



SS 



+ 



W V 

SS SS 



SS 



— (2At;2,m + 0 






>2 



SS 



SS 






w„ a/ 1 + r 



SS 



w + 



+ 1) 



SS 



SS 



Vi 



w„ a/ 1 + r 



,/2 



(w-)' 



SS 



+ 



1-2A- 






w 



SS 



^22 



A 



. ViTtf _ 






(4-19) 



The subscript ss indicates the steady state solution and the superscript * 
denotes the perturbation. These equations are for the case of fixed A and Vf. Thus, 
the corresponding conditions, equation (3-21) and v =0 at z = f ^pply to this case. 

Furthermore, the axial force at the freeze line T^ is allowed to fluctuate whose 
linearized equation is given as 



2 ( 1-0 



SS 



Qw* = [ 2(1 - f){w„v;, -I- w;v„ } -I- 2fw„<7, , „ ]r* 



+[2(1 - f)rX + 2fr„(T„ „ ]w’ -t- [2(1 - f)r„v„ ](w’) 
+[2(l-f)r„w;,]v*+[2(l-f)r„w„](v*) +[2fr„w„](r;, - T; 



(4-20) 



74 



In addition to the typical boundary conditions for r*,w* and v* at the origin and at the 

freeze line which are given in equations (3-18) and (3-19), the following initial stress 
conditions apply for the two-layer flow: 




at z = 0 



(4-21) 



4.4. Computational Procedure 



The computational procedure for the two-layer case is in principle the same as 
that of Newtonian single-layer case although its formulation is much more 
complicated and it requires more computation time. The linearized equations (4-14) 
through (4-20) are converted into the following algebraic equations using a finite 
differencing scheme. 




(4-22) 




(4-23) 






(4-24) 



75 



/' ■> 

* 

r 




* 

r 


3= 


w‘ 


1 ^ 


^ = C4^ 


> 




Sn 


*11 




T # 


V. 




^22. 






f 4 {v'} 



( 4 - 25 ) 






r 




r* 


= 


w* 


i ^ 


► = 65^ 


> 


— * 




S33 


^33 






V. V 




r 






f 5 {v-} 



( 4 - 26 ) 








r • 1 

r 


fw’ 


« 


- r = 


14 J 


t* 

[* 22 j 








f 6 {v'} 



( 4 - 27 ) 



Here the components of the vectors r*, w',s*,,S3j and represent the (n- 1 ) unknown 

r 1 T X 

values of each variables at the nodal points (i.e., r’ =[r2,---,r„' ,s,‘, =[^n2>‘“>^nn] > 

*33 =[<733,2 >■ ■•>*^33,11] .*11 =[^ii,2>‘".^n4i] . etc.) where n corresponds to the freeze line 

height. The first three equations are from the continuity equation, the axial and the 
radial force balances, respectively and the last three are from the constitutive relations 
for the normal stress components. By combining the continuity equation, the axial 
force balance and the three constitutive relations, the following equation is obtained. 



DiQx' = Eix’ + Fiv* 



( 4 - 28 ) 



* [’’2 '''^2 >’".'''d.< 7 i|, 2 .‘".< 7 ,,„,C 733 2 ,”',CT 33 j,,T 22 _ 2 ,’”, ^22 n ] 



and 



v’ =[v 2 .v;,---,v;_,,t;]'^ 



where 



76 



Here \ and v* are the eigenvectors of (5n-5) and (n-1) dimensions, respectively. The 
matrices Di and Ei are (5n-5)x(5n-5) matrices and Fi is a (5n-5)x(n-l) matrix. The 
radial force balance (4-24) along with the expression for T^’ (equation (4-20)) result in 

DzQx* = E2 x’ + F 2 V* f4-29) 

Here D 2 and E 2 are (n-I)x(5n-5) matrices and F 2 is a (n-l)x(n-l) matrix. By 
multiplying the inverse of F 2 to equation (4-29), an explicit expression for v' is 
obtained. By substituting this expression for v* into equation (4-28) the standard form 
for an algebraic eigenvalue problem is obtained which can be solved using the 
subroutine GVLRG in IMSL package.^^ 

Since the eigenvector x* has (5n-5) unknown components, excessive time is 
required for the eigenvalue problem even for a relatively low spatial resolution. When 
n=41 (i.e., 41 data points between the origin and the freeze line height), the matrix 
size that should be inverted is 200x200. The dependence of the first eigenvalue with 

the largest real part on the number of discretization is shown in Table 4-1. 
Eigenvalues for n— is obtained by the use of Richardson extrapolation method. 
Two extrapolation values for different numbers of data points are given in the table. 
The first one is computed with the first three data points (n=21, 31 and 41) whereas 
the second one is with all the five data points (n=21,31,41,51 and 61). Although 
there is a slight difference between the two extrapolation values, most of the results 
presented in this chapter are based on the extrapolation with the three data points in 
order to avoid the excessive computation time. 



77 



Table 4-1 Dependence of the first eigenvalue on the number 

of discretization (The values for n->®« have been 
calculated by the Richardson extrapolation method; 

i=5, f=0.2, A=0.1, T,=3.23, B=0.2, BUR=3.63 
and tr=54.35) 



n 


Qr 


Cli 


21 


-1.0823 


18.89 


31 


-0.7478 


30.08 


41 


-0.6552 


41.25 


51 


-0.6050 


52.34 


61 


-0.5772 


63.38 


n— with the 
three data point 


-0.5547 


60.32 


n— with five 
data point 


-0.5724 


146.12 



78 



4.5. Results and Discussion 



The results of the stability analysis for the two-layer film are plotted on BUR 
vs. tf plane as in the previous chapter. 

In Figure 4-1, the neutral stability curve is described on the constant B 
contours of the steady state solution. The outer-layer flow rate fraction f is 0.2, the 
dimensionless relaxation time is 0.1, and the freeze line height is 5. Also specified in 
the figure are the real and the imaginary part of the eigenvalues for several data 
points. The constant B contours for this f=0.2 case are not much different from those 
of a Newtonian single-layer case, (see Figure 2-7 or Figure 3-5 for comparison) The 
overall stability behavior, however, is different in that the upper unstable region (i.e., 
the unstable region where BUR>3.5) of the Newtonian single-layer flow (Figure 3-5) 
is not obtained in the two-layer coextrusion flow. 

In Figure 4-2, the constant B contours and the neutral stability curve for a 
larger value of f (i.e., f=0.8) are given. The constant B contours for this case 
resemble those of a UCM single-layer flow, indicating a strong influence of the 
viscoelastic outer-layer (see Figure 2-7). The neutral stability curve is shifted to the 
left compared to that of Figure 4-1 for f=0.2. Thus, the region of stable operation is 
decreased as the thickness of the viscoelastic layer is increased. 

The neutral stability curves for various values of f are collected in Figure 4.3. 
Also given in the figure are the neutral stability curves for the single-layer flows of a 
Newtonian (f=0.0) and an UCM fluid (f=1.0) for comparison. The dimensionless 
freeze line height is fixed at five, and the dimensionless relaxation time is 0. 1 for all 
the data. The most interesting aspect of these neutral stability curves is that the 
Newtonian single-layer film is unstable when BUR is larger than about 3.5 whereas 
the upper unstable region is not observed for all other cases. For all cases there exist 



79 




Thickness Reduction 



Figure 4-1 The first eigenvalue of the two-layer case at fixed A and 

Vf(The first and the second number in the parenthesis 
represent the real and the imaginary part of the 
eigenvalue; Dashed line is a neutral stability curve at the 
left hand side of which the bubble is unstable; 
dimensionless freeze line height ^ =5; dimensionless 

relaxation time X=0. 1; flow rate ratio=0.2) 



80 




Thickness Reduction 



Figure 4-2 The first eigenvalue of the two-layer case at fixed A and 

Vf(The first and the second number in the parenthesis 
represent the real and the imaginary part of the 
eigenvalue; Dashed line is a neutral stability curve at the 
left hand side of which the bubble is unstable; 
dimensionless freeze line height f =5; dimensionless 

relaxation time X=0. 1; flow rate ratio=0.8) 



81 




Thickness Reduction 



Figure 4-3 The effect of flow rate ratio f on the stability of the two- 

layer film at fixed A and Vf (The curves are neutral stability 
curves; dimensionless freeze line height f =5; 

dimensionless relaxation time X=0. 1) 



82 




Figure 4-4 The effect of flow rate ratio f on the stability of the two- 

layer film at fixed A and Vf (The curves are neutral stability 
curves; dimensionless freeze line height £ =5; 

dimensionless relaxation time X=0.05) 



83 




Figure 4-5 The effect of flow rate ratio f on the stability of the two- 

layer film at fixed A and Vf (The curves are neutral stability 
curves; dimensionless freeze line height t =5; 

dimensionless relaxation time X=0.2) 



84 



a maximum film thickness reduction beyond which the bubble becomes unstable. The 
maximum film thickness reduction appears to be increasing with the blow-up ratio 
(BUR). As the outer-layer flow rate fraction is decreased, the neutral stability shifts 
to the right, broadening the region of stable operation. Thus, it is indicated that 
smaller fraction of viscoelastic layer is more favorable for the stability of the bubble. 

In Figure 4-4 and 4-5, similar neutral stability curves are given for two 
different values of the dimensionless relaxation time of the outer-layer X. Compared 
to Figure 4-3, the general trend remains the same. For the smaller value of X, 
however, the neutral stability curves for f>0.2 shift to the right and all are located on 

the right hand side of the Nev/tonian curve. When X is large (i.e., X=0.2), on the 
other hand, the curves shift to the left and all are now located on the left hand side of 
the Newtonian curve. These trends may be conceivable considering the fact that the 
viscoelasticity tends to limit the elongational deformation of a material. Thus, the 
large relaxation time X acts against the stretching of the film, resulting in a decrease 
in the maximum thickness reduction for a stable operation. 

4.6. Summary and Conclusion 

The stability of a two-layer blown film coextrusion has been investigated by 
the method of a linear stability analysis. Only the axisymmetric disturbances have 
been considered for fixed A (amount of air trapped in the bubble) and fixed Vj- (the 

take-up velocity at the freeze line). The overall stability behavior may be summarized 
as follows: 

(1) The presence of the viscoelastic layer removes the upper unstable region of the 
Newtonian single-layer flow, thus making a stable operation possible at a higher value 



85 



of BUR. This stabilizing influence in terms of maximum BUR is apparent regardless 
of the relaxation time and the outer-layer flow rate fraction f 

(2) The maximum film thickness reduction for a stable operation is predicted to 
increase with BUR for all cases. 

(3) The location of the neutral stability curves varies depending on the magnitude of 
the relaxation time. As X decreases, the neutral stability curves shift to the right and 
broaden the region of stable operation. 

(4) For a fixed value of X, the neutral stability curves moves to the right with 
decreasing f (the flow rate fraction of the viscoelastic outer-layer). Thus, a smaller 
value off is always favorable for the bubble stability. 

(5) When X is moderately small (e.g., X=0.05), the neutral stability curves for the 
coextrusion flow (0.2 <f <0.8) are located on the right hand side of either a 
Newtonian (f=0.0) or a UCM single-layer (f=1.0) flows. Thus, a kind of synergism in 
terms of bubble stability is predicted with the coextrusion flow. For example, when 
A, =0.05 and f=0.2, a stable operation is possible at BUR=3.0 and t^=80. In case of 

single-layer flows of either a Newtonian or a UCM fluid, however, such a high 
thickness reduction at BUR=3.0 is not possible. 

Even with the simple model presented in this thesis, a very complicated 
stability behavior is predicted for a blown film process. Although a thorough 
physical explanation for such behaviors can not be given at present, many of the 
predictions appear to be in reasonable agreement with experimental observations 
at least qualitatively. It is interesting to point out that the upper unstable region 
was in fact observed with a LLDPE although the instability mode in experiment 
was a helical instability rather than an axisymmetric one. The stabilizing influence 
of a viscoelastic layer (especially the increase in the maximum obtainable BUR) is 



86 



also observed frequently in LLDPE/LDPE blown film coextrusion. Due to the 
simplicity of the present model, quantitative agreement between experiment and 
theory can not be expected. The present study, nevertheless, may serve as the 
seminal work for future studies to apply more complicated but realistic models to 
this extremely complicate process. 



CHAPTER 5 
EXPERIMENT 

5. 1. Introduction 

Experiment on a bicomponent two-layer blown film coextrusion has been 
conducted using a linear low density polyethylene (LLDPE) and a low density 
polyethylene (LDPE) as the inner and the outer layer, respectively. These two 
materials are chemically the same, but have distinctly different rheological 
properties due to the differences in their molecular structure. Furthermore, the 
physical properties of the films manufactured using these materials are also very 
different from one another due to a similar reason. The major objective of the 
experiment is to investigate the physical properties of the film as influenced by the 
coextrusion process. Films with various layer structures (i.e., different values of 
f) have been fabricated and their tensile and tear properties have been measured 
and compared. During the extrusion, the processing characteristics (especially the 
bubble stability) of the two-layer films have been also investigated. While the 
flow characteristics of the extrusion process may be studied through a process 
modeling as presented in the previous chapter, the study on film physical 
properties is strictly empirical at the current stage of knowledge. Since the film 
physical properties and the processing characteristics of the composite films are 
important factors for determining the choice of materials for the coextrusion 
process, both the process modeling and the empirical evaluation of the film 



87 



88 



physical properties should be the integral part of the study on the blown film 
coextrusion process. 

The results indicate that the addition of a LDPE layer to a pure LLDPE 
layer improves the bubble stability as expected. However, the tensile properties 
of the composite film are shown to have a detrimental influence by the 
coextrusion. Although the present experimental study is limited in its scope, the 
results suggest that a precaution should be taken in selecting the materials for 
coextrusion so that a critical deterioration of important physical properties can be 
prevented. 



5.2. Experimental 



5.2.1. Materials 

A low density polyethylene (LDPE) and a linear low density polyethylene 
(LLDPE) have been used as the inner and the outer layer, respectively. Both 
materials are commercial products of Union Carbide Corporation whose product 
names are DYNK-9 and GRSN-7047, respectively. 

The LDPE (i.e., DYNK-9) is produced by the conventional high pressure 
process whereas the LLDPE (GRSN-7047) is produced by a new low pressure 
process. Both materials are of 1 melt index and 0.918 density. Due to the 
differences in the manufacturing process, the two materials have distinctly 
different molecular structure in that the LDPE has a high degree of long chain 
branching whereas the LLDPE does not. Consequently, the LDPE has a larger 
relaxation time than the LLDPE and shows a pronounced strain hardening 
behavior which is usually absent with the LLDPE. Due to these rheological 

W 

differences, the LDPE has much better bubble stability although the physical 
properties of its films are much inferior to the LLDPE. 



First Normal Stress (N/m 



89 



IN 




Figure 5-1 First normal stress difference (x,i-T 22 ) of LDPE and LLDPE 

evaluated using plate-plate viscometer (data for LDPE were 
measured at 145°C and data for LLDPE at 200°C) 



Shear Viscosity (Kg/m-sec 



90 




Figure 5-2 Shear viscosity of LDPE and LLDPE evaluated using plate- 

plate viscometer (° and □) when y < Is”' and capillary 

viscometer (• and ■) when y > Is”' (data for LDPE were 
measured at I45°C and data for LLDPE at 200°C) 



91 



In Figure 5-1 and Figure 5-2, the first normal stress difference (t,,-T 22 J and 
the shear viscosity at different shear rates are plotted. In Figure 5-1, the viscosity 
measuring temperature (145°C) is different from the experiment temperature 
(160°C) for LDPE. The temperature dependence of the first normal stress 
difference, however, is reported to be relatively small. Compared with the first 
normal stress difference measured at 145°C, that at 160°C decreases by less than 
5%. Thus, although the temperature dependence is considered, Figure 5-1 shows 
that the first normal stress difference of LDPE is larger than that of LLDPE at the 
strain rates which are within the range of experiment. It indicates that LDPE is 
more viscoelastic than LLDPE. In Figure 5-2, it is observed that both LDPE and 
LLDPE have similar viscosity in the shear rate of the experiment. 

5.2.2. Extrusion Conditions 

Two 19 mm single screw extruders were used separately for the inner and 
the outer layer. The diameter and the gap of the coextrusion die were 25.4 mm 
and 1 mm, respectively, and a single lip air ring was used for bubble cooling. It 
was attempted to make films with an average thickness of about 30 pm for five 

different flow rate ratios; LLDPE/LDPE=100/0,80/20, 50/50, 20/80, and 0/100. 
The details of the operating conditions are summarized in Table 5-1. 

5.3. Results and Discussion 
5.3.1. Processing Characteristics 

The differences in bubble shapes were difficult to detect since the bubble 
diameters were very small. The bubble stability between the samples, however, 
could be compared in terms of maximum obtainable blow up-ratio (BUR). The 



92 



Table 5-1 Processing conditions 



Sample No. 


1 


2 


3 


4 


5 


LLDPE/LDPE 


10/0 


8/2 


5/5 


2/8 


0/10 


Extruder 1 for the LLDPE ('inner-laver') 




RPM 


55 


56 


28 


10 


0 


Extruder 2 for the LDPE fouter-laver) 




RPM 


0 


12 


30 


48 


60 



Total output 



(g/min) 


19 


26 


20 


27 


27 


Blow up ratio 


1.7 


1.7 


1.9 


1.9 


1.9 


Film gauge (pm) 


24 


39 


28 


27 


36 



Melt temperature in Extruder 1= 200 °C 
Melt temperature in Extruder 2= 160 °C 
Die temperature= 220°C 



93 



BUR is one of the most important parameters which may have the strongest 
influence on the physical properties of the film. 

It was attempted throughout the experiment to maintain the BUR at about 
2.0. For samples 3, 4 and 5 which consist of more 50% of the LDPE, the BUR of 
1.9 could be easily obtained. For samples 1 and 2 for which the LLDPE is the 
major constituent, however, the bubbles were much more unstable and the highest 
BUR that could be obtained was only 1.7. In case of sample 1 (LLDPE single- 
layer bubble), even the extrusion output had to be lowered to obtain a stable 
bubble at the BUR of 1.7. This result indicates that the poor bubble stability of 
LLDPE can be improved by coextruding it with a more stable LDPE. These 
observations appear to be a qualitative agreement with the theoretical prediction 
described in the previous chapter. 

5.3.2 Physical Properties 

The physical properties of the film samples that have been evaluated 
include the tensile properties (tensile strength, % elongation and yield strength) 
and the tear resistance (Elmendorf tear). These properties are the most important 
ones that are routinely measured in industry for the evaluation of the film physical 
properties, and represent two different types of tests. 

The tensile strength measurement is a low strain rate test in which the 
response of the film samples to a slow deformation is measured. In the tear 
strength measurement, on the other hand, the rate of deformation is very fast. 
Both test methods are standard ASTM methods (tensile strength measurement; 
ASTM D412-51T; tear strength measurement; ASTM D689-44) and the 
measurements were done by the physical property testing laboratory of Union 
Carbide Corporation. In Figure 5-3, the shape of the specimen and the direction 
of the deformation are described schematically for both tests. 



94 



~ 1 .25" 




Direction of deformation (rate; 20" /min) 



(a) 



2.5" 




Direction of deformation 
(perpendicular to sample) 



(b) 



Figure 5-3 Specimens for the measurement of physical properties 

(a) specimen for tensile property measurement 

(b) specimen for Elmendorf tear measurement 



95 



In the tensile property test, the force per unit area (e. g., MPa) that is 
required for the given rate of deformation is measured as a function of strain until 
the failure (or breakage) of the sample occurs. From this continuous 
measurement, the yield strength (for a ductile deformation), tensile strength (at 
break) and the % elongation at break are determined. The tear resistance is 
represented by the force per unit thickness of the film (e.g., g/mil) acting against 
the fast deformation. The direction of the applied force and the direction of the 

tear propagation are perpendicular to each other in this test as described in Figure 
5-3. 

In both measurements, test specimens in the machine direction (or axial) 
and the transverse (or circumferential) direction are taken for measurements. The 
samples in the machine direction (MD) and in the transverse direction (TD) show 
different values since the deformation history of the film in the two directions are 
different. Since the deformation in the machine direction is more severe during 
the film blowing process, the tensile strength in the machine direction is typically 
higher than that in the transverse direction. The imbalance of the tear resistance 
in the machine and the transverse direction represents the anisotropy in the film 
physical properties which is influenced by both processing conditions and the 
intrinsic properties of the material being processed. The results of the physical 
property measurements are summarized in Table 5-2 and described graphically in 
Figures 5-4, 5-5 and 5-6. The physical properties of the single-layer LLDPE and 
the LDPE films are well known and the results for the film samples 1 and 5 are in 
good agreement with the known values. It is noteworthy that the machine- 
direction (MD) tensile strength of the composite films (samples 2,3, and 4) is as 
low as that of the LDPE film. It seems that the "mutual interlayer destruction” is 
occurring for these samples. 26 As can be seen from Table 5-2 and Figure 5-4, the 
tensile strength of the single-layer LDPE (sample 5) is much weaker than that of 



96 



Table 5-2 Tensile and tear properties of the samples 



Sample number 1 2 3 4 



LLDPE/LDPE 


10/0 


8/2 


5/5 


2/8 


Blow up ratio 


1.7 


1.7 


1.9 


1.9 


Gauge (lim) 
Tensile strength 
(MPa) 


24 


39 


28 


27 


MD 


38 


25 


24 


26 


TD 

Yield strength 
(MPa) 


22 


21 


21 


18 


MD 


14 


- 


- 


- 


TD 

Elongation at 

break (%) 


9 


9 


10 


10 


MD 


550 


240 


180 


120 


TD 

Elmendorf tear 
(g/mil) 


660 


680 


670 


650 


MD 


33 


24 


365 


290 


TD 


260 


127 


66 


74 



_5 

0/10 

1.9 

36 

28 

20 



11 

160 

500 

101 

73 



(lmil= 0.001 inch) 



Tensile strength (MPa) 



97 




10/0(1) 8/2(2) 5/5(3) 2/8(4) 0/10(5) 

LLDPE/LDPE ratio (sample number) 



m MD MTD 






Figure 5-4 Tensile strength of the samples (The number in the parenthesis 

is the sample number) 




Elongation at break (%) 



98 



800 

700 

600 

500 

400 

300 

200 

100 

0 




10/0(1) 8/2(2) 5/5(3) 2/8(4) 0/10(5) 

LLDPE/LDPE ratio (sample number) 




Figure 5-5 Elongation of the samples (The number in the parenthesis 

is the sample number) 



Elmendorf tear (g/mil) 



99 




10/0(1) 8/2(2) 5/5(3) 2/8(4) 0/10(5) 

LLDPE/LDPE ratio (sample number) 




MD MID 



Figure 5-6 Elmendorf tear of the samples (The number in the parenthesis 

is the sample number) 



100 



LLDPE (sample 1). Furthermore, it is less ductile than LLDPE as it has no yield 
strength and low elongation (Figure 5-5). Consequently, fractures which are 
formed first in the LDPE layer may propagate into the LLDPE layer, resulting in a 
low MD tensile strength for the composite films. Since the two materials are 
chemically compatible, a good adhesion between the two layers is virtually 
ensured and it may not be surprising to see such trend in tensile strength (i.e., 
mutual interlayer destruction). It seems that no yield strength in machine 
direction and the low elongation at break of the composite films also support the 
argument. The tear resistance, on the other hand, may be less affected by the low 
ductility of the LDPE since it is a very high deformation rate testing. Although 
the results of the measurements are rather erratic, it is interesting to note that the 
MD tear strength of the composite films (samples 3 and 4) shows much higher 
values than those of the single-layer films of either LLDPE or LDPE. This result 
may suggest the "mutual enforcement" of the physical property in this high strain 
rate test.^6 

5.4. Summary and Conclusion 

In the present experiment in which an LLDPE and an LDPE were used, 
enhancement in bubble stability was observed with the composite structures 
compared to that of the single-layer LLDPE. In case of the physical properties of 
the film, however, a significant deterioration of the MD tensile strength was 
observed, whereas a drastic increase was observed in the MD tear resistance. It 
may be due to the combined influence of the low ductility and the low toughness 
of the LDPE that results in the propagation of the craze formed in the LDPE layer 
into the LLDPE layer at a relatively low tensile stress. The resistance of the 



101 



composite film to a high strain rate deformation, however, appear to be dictated 
by a different mechanism in that a mutual enhancement is observed. 

While the coextrusion process is used to provide improved end-use 
properties to the products through collective improvement, this study suggests 
that some properties can be significantly deteriorated by coextrusion unless the 
materials to be coextruded are chosen carefully. 



CHAPTER 6 

SUMMARY AND CONCLUSION 



In the present study, a bicomponent two-layer blown film coextrusion has been 
investigated both theoretically and experimentally. The objective of the theoretical 
analysis is to investigate the competing influence of two rheologically different fluids on 
the bubble dynamics by the use of a simple mathematical model. An isothermal two-layer 
flow is assumed, in which the inner and the outer layers consist of a Newtonian fluid and 
an Upper-convected Maxwell fluid, respectively. For the derivation of governing 
equations, "thin film approximation" was applied, in which the velocity and the stress 
fields are assumed to be uniform across the film thickness. 

In the steady state analysis, the governing equations are solved numerically to 
determine the bubble radius and the film thickness profile as a function of material 
properties and processing conditions. When the relaxation time is small, the bubble 
dynamics is not much different from that of a Newtonian single-layer flow. Due to the 
existence of the Maxwell fluid layer, however, it is found that with increasing relaxation 
time, the viscoelastic effect becomes more pronounced and eventually dominates the two- 
layer bubble dynamics even at a smaller flow rate fraction of the Maxwell fluid layer. 

In the stability analysis, both Newtonian single-layer and the two-layer cases are 
investigated. An operating condition of a fixed amount of air (A) and take-up velocity (Vf) 
are mainly considered for its relevance to the industrial practice. By the use of linear 
stability analysis, neutral stability curves have been determined on a BUR-t^ plane for the 
range of practical interest(i.e., l^UR^ and t,^0). In case of the Newtonian single- 
layer flow, the computational results predict that the bubble stability is limited by both 



102 



103 



BUR and regardless of the freeze line height, which is in quantitative agreement with 
experimental observations for linear low-density polyethylene's. It also indicates that there 
exists an optimum freeze line height for a maximum stable operating region. In case of 
two-layer flow, on the other hand, the upper unstable operating region appearing in the 
Newtonian single-layer case for large BUR, does not exist, indicating the stabilizing 
influence of the viscoelastic layer in terms of BUR. With respect to the thickness 
reduction, however, the stable operating region can be larger or smaller depending on the 
relaxation time as well as the flow rate fraction of the Maxwell fluid layer. Thus, it can be 
concluded that there exist an optimum relaxation time and an optimum flow rate fraction 
of the Maxwell fluid layer for a maximum stable operating region. 

LLDPE and LDPE were used as the inner and the outer layers in the experimental 
study. It is known that LDPE has a larger relaxation time than LLDPE, thus having better 
stability than LLDPE in terms of BUR. As expected, LLDPE bubbles were stabilized by 
the coextrusion with an LDPE, indicating the stabilizing influence of LDPE. It was also 
observed that the tensile properties of LLDPE were deteriorated by coextrusion while the 
tear strength was improved. These trends in film physical properties are not predictable 
theoretically although they are related to processing conditions to some extent. Thus, it 
may be concluded that the bubble stability of LLDPE can be improved by coextrusion with 
materials of an optimum relaxation time, but the materials should be chosen carefully to 
prevent a deterioration of the physical property of the coextruded film. 



APPENDIX A 
COMPUTER PROGRAM 
FOR STEADY STATE SOLUTION 



The computer programs for steady state solution of the two layer film and the 
Maxwell fluid single layer film are presented. The program for the Newtonian single 
layer film is not listed here since it has no difference from the works that have been 
done before. The fourth-order Runge-Kutta method is used for integration, and the 
operating parameters and the boundary conditions used are listed in Chapter 2. The 
subprogram FUNCTION RUNGE for integration is listed at the end of the main 
program for the two-layer film, and omitted in the Maxwell fluid single layer case. 
The direction of integration is from the die exit to the freeze line. 

A-A. 1 Computer Program For Steady State Solution 

Of The Two-Laver Film 



P ****^4*:^*t*****^l^^*i^t ****************************************** 

C COMPUTER PROGRAM FOR STEADY STATE SOLUTION OF 
C THE TWO-LAYER BLOWN FILM COEXTRUSION 

C BY FOURTH-ORDER RUNGE-KUTTA METHOD 

C WITH DOUBLE PRECISION 

p ************************************************************** 

C 

IMPLICIT REAL*8 (A-H.O-Z) 

INTEGER RUNGE 
DIMENSION F(6),R(6) 

DATA Z,DZ,ICOUNT,IFREQ/O.DO, 1 .D-3,0,200/ 

C 

C * ♦ * * DESCRIPTION OF VARIABLES * • * * 

C 

C R(l)=r R(2)=r' R(3)=w 

C F(l)=r’ F(2)=r" F(3)=w’ 

C R(4)=0|[ R(5)=Oj 3 R(6)=t22 

C F(4)=o,i' F(5)=033- F(6)=l22‘ 

C 

C •♦‘♦PARAMETERS**** 



104 



105 



C 

C T= T-BAR, B= INFLATING PRESSURE, 

C FR= FLOW RATE RATIO, FL= FREEZE LINE HEIGHT, 

C VISCO= VISCOSITY RATIO, D= RELAXATION TIME. 

C 

C **** INITIAL GUESS OF R-PRIME ♦*** 

C 

PRINT *,'INPUT D,FR,B,T & RPI’ 

READ (*,*)D,FR,B,T,RPI 

VISCO=I.DO 

FL=5.D0 

G=(I.DO-FR)*VISCO/FR 
WRITE (*,104) T,B,D,FR,VISCO,FL 
IP=0 
C 

C ♦♦*• INITIAL VALUES **♦* 

C 

3 Z=0.D0 
R(I)=1.D0 
R(2)=RPI 
R(3)=1.D0 
R(4)=0.ID0 
R(5)=0.ID0 
R(6)=0.D0 
C 

1 IF (IP.EQ.O) GO TO 2 

WRITE(*,I02) Z,R(1),R(2),R(3),R(4),R(5),R(6) 

C 

C ♦♦** FOURTH-ORDER RUNGE-KUTTA METHOD**** 

C 

2 K=RUNGE(6,R,F,Z,DZ) 

IF(K.NE.l)GOTO 10 

TI=R(1)*R(3)*DSQRT(I.D0+R(2)*R(2))/D 

T2=R(2)/R(I) 

T3I=-T2/2.D0+TI*D*R(4)/(2.D0*G) 

T3=T3I-(I.D0+R(2)*R(2))*(T+B*R(I)*R(I))/(4.D0*FR*G) 

F(I)=R(2) 

F2I=(R(5)+G*(T2-T3)/(TI*D))/(R(4)-G*(T2+2.D0*T3)/(TI*D)) 

F(2)=(1.D0+R(2)*R(2))*(-2.D0*B*R(1)/(T+B*R(I)*R(I))+F2I/R(I)) 

F(3)=T3*R(3) 

F4I=-(TI+2.D0*(T2+T3))*R(4)-2.D0*(2.D0*T3+T2)*R(6) 

F(4)=F4I-(2.D0*T3+T2)/D 

F(5)=(2.D0*T2-T1)*R(5)+(T2-T3)*(2.D0*R(6)+I.D0/D) 

F(6)=T3/D+(2.D0*T3-T1)*R(6) 

GO TO 2 
C 

C **** PRINT INTERVAL **** 

C 

10 IF (Z.GE.FL) GO TO 4 
ICOUNT= ICOUNT+I 
IF (ICOUNT.NE.IFREQ) GO TO 2 
ICOUNT=0 
GO TO I 
C 



106 



C **** TERMINATION 
C 

4 RPF=R(2) 

IF (DABS(RPF).LE.l.D-3) GO TO 5 
WRITE (*,103) RPI,RPF 
RPI=RPI+5.D-2 
GO TO 3 

5 IF (IP.NE.l)GOTO 222 
TZ=T+B*R(1)**2 
TR=1.D0/R(3) 

WRITE (*,221) TZ,B,R(1),TR 

221 FORMAT(/13H DRAW FORCE= ,F9.3,5X,4H B= ,F6.3/ 

&6H BUR= ,F6.3,12X,5H Tr= ,F9.3/) 

STOP 
222 1P=1 

WRITE (*,104) T,B,D,FR,VISCO,FL 
WRITE (*,101) DZ 
GO TO 3 

101 FORMAT (/5H DZ= ,F7.5/2X,3H X ,8X,1HR,7X,7HR-PR1ME, 
&6X,1HW,8X,5HSIG11,6X,5HSIG33,6X.5HTAU22) 

102 FORMAT (F7.3,5(1X,D10.3),1X,D9.3) 

103 FORMAT (6H RPI=,D13.5,6H RPF=,D13.5) 

104 FORMAT(/7H T-BAR=,F9.3,5H B=,F6.3,9H LAMDA=,F6.3/ 
9HF-RATIO=,F6.3,10H V-RATIO=,F6.3,13H FR0ST-L1NE=,F6.3/) 

END 

C 

C 

C 

C ******** SUBPROGRAM FUNCTION RUNGE ********* 

C (FOURTH-ORDER RUNGE-KUTTA METHOD) 

C 

C 

FUNCTION RUNGE(N,Y,F,X,H) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

REAL*8 Y,F,X,H 
INTEGER RUNGE 

DIMENSION PHI(50),SAVEY(50),Y(N),F(N) 

DATA M/0/ 

M=M+1 

GOTO(l,2,3,4,5),M 

C 

1 RUNGE=1 
RETURN 

C 

2 DO 22 J=1,N 
SAVEY(J)=Y(J) 

PHI(J)=F(J) 

22 Y(J)=SAVEY(J)+.5D0*H*F(J) 

X=X+.5D0*H 

RUNGE=I 

RETURN 

: 

3 DO 33 J=1,N 
PHI(J)=PHI(J)+2.D0*F(J) 



non 



107 



33 Y(J)=SAVEY(J)+.5DO*H*F(J) 
RUNGE=1 
RETURN 

C 

4 DO 44 J=1,N 
PHI(J)=PHI(J)+2.D0*F(J) 

44 Y(J)=SAVEY(J)+H*F(J) 

X=X+.5D0*H 

RUNGE=1 

RETURN 

5 DO 55 J=1,N 

55 Y(J)=SAVEY(J)+(PHI(J)+F(J))*H/6.D0 
M=0 

RUNGE=0 

RETURN 

END 



A-B.2 Computer Program For Steady State Solution 
Of Maxwell Fluid Single-Laver Film 



C COMPUTER PROGRAM FOR STEADY STATE SOLUTION OF 
C BLOWN FILM EXTRUSION OF MAXWELL FLUID SINGLE-LAYER 
C BY FOURTH-ORDER RUNGE-KUTTA METHOD 
C WITH DOUBLE PRECISION 

IMPLICIT REAL*8 (A-H.O-Z) 

INTEGER RUNGE 
DIMENSION F(6),R(6) 

DATA Z,DZ,ICOUNT,IFREQ/0.D0,l.D-3, 0,200/ 

C 

C R(I)=r R(2)=r' R(3)=w 

C F(l)=r' F(2)=r* F(3)=w' 

C R(4)=an R(5)=033 R(6)=i22 

C F(4)=aii' F(5)=a33' F(6)=T22* 

C 

C PARAMETERS 
C 

C T= T-BAR, B= INFLATIO PRESSURE, 

C D= RELAXATION TIME, FL=FREEZE LINE HEIGHT 
C 

D=0.ID0 

FL=5.D0 



♦♦♦♦ INITIAL GUESS FOR R-PRIME 

READ (♦,♦) B,T,RPI 
WRITE (*,104) T,B,D,FL 
IP=0 



^ ^ ^ CJCJCJU CJCJCJ UCJCJ 



108 



INITIAL VALUES **** 



3 Z=0.D0 
R(I)=I.D0 
R(2)=RPI 
R(3)=I.D00 

R(4)=DSQRT(1.+R(2)**2)*(T+B*R(I)**2)/(2.*R(I)*R(3)) 

WPR=-(I.+RPI**2)/4.*(T+B)-RPI/2. 

EE=I./DSQRT(1.+R(2)**2) 

RNEI=-(2.*WPR+RPI)*EE 

RNE2=WPR*EE 

RNE3=(RPI-WPR)*EE 

R(5)=RNE3 

R(6)=RNE2 

1 IF (IP.EQ.O) GO TO 2 

WRITE(6,102) Z,R(1),R(2),R(3),R(4),R(5),R(6) 

•♦♦♦ FOURTH-ORDER RUNGE-KUTTA METHOD ***♦ 

2 K=RUNGE(6,R,F,Z,DZ) 

IF (K.NE.l)GOTO 10 
AA=DSQRT(1.+R(2)**2) 

F(1)=R(2) 

F(2)=AA**2*(-2.*B*R(1)/(T+B*R(1)**2)+R(5)/(R(1)*R(4))) 

RD=F(2) 

F{3)=R(3)/(R(4)+4.*R(6)+2./D)*(-R(2)*RD* 

&(T+B*R(1)**2)/(2.*R(1)*R(3)*AA)-B*R(2)* 

&AA/R(3)-(R( 1 )*R(3)* AA/D+R(2)/R( 1 ))*R(4) 

&- 2 . *R(2)*R(6)/R( 1 )-R(2)/(R( 1 )*D)) 

WP=F(3) 

BB=WP/R(3) 

CC=R(2)/R(1) 

F(4)=-(R(1)*R(3)*AA/D+2.*(BB+CC))*R(4) 

&-2.*(2.*BB+CC)*R(6)-(2.*BB+CC)/D 

F(5)=(2.*CC-R(1)*R(3)*AA/D)*R(5)+ 

&2.*(CC-BB)*R(6)+(CC-BB)/D 

F(6)=BB/D+(2.'*BB-R(1)*R(3)*AA/D)*R(6) 

GO TO 2 

**♦* PRINT INTERVAL *♦♦♦ 

10 IF (Z.GE.FL) GO TO 4 
1C0UNT= ICOUNT+ 1 
IF (ICOUNT.NE.IFREQ) GO TO 2 
ICOUNT=0 
GO TO 1 



**♦* TERMINATION 



4 RPF=R(2) 

IF (DABS(RPF).LE.I.D-2) GO TO 5 
WRITE (*,103) RPI.RPF 



109 



RPI=RPI+5.D-2 
GO TO 3 

5 IF (IP.NE.l)GO TO 222 
TZ=T+B*R(1)**2 
TR=1.D0/R(3) 

WRITE (6,221) TZ,B,R(I),TR 

221 FORMAT(/13H DRAW FORCE= ,F9.3,5X,4H B= ,F6.3/ 

&6H BUR= ,F6.3,12X,5H Tr= ,F9.3/) 

STOP 
222 IP=1 

WRITE (6,104) T,B,D,FL 
WRITE (6,101) DZ 
GO TO 3 

101 FORMAT (/5H DZ= ,F7.5/2X,3H X ,8X,1HR,7X,7HR-PR1ME, 

6 6X,1HW,8X,5HSIG1 1,6X,5HSIG33,6X,5HTAU22) 

102 FORMAT (F7.3,5(1X,D10.3),1X,D9.3) 

103 FORMAT (6H RPI=,D13.5,6H RPF=,D13.5) 

104 FORMAT(/7H T-BAR=,F9.3,5H B=,F6.3,9H LAMDA=,F6.3/ 
& 13H FROST-LINE=,F6.3/) 

END 



APPENDIX B 

COMPUTER PROGRAM FOR STABILITY ANALYSIS 

The computer programs for stability analysis of Newtonian single-layer film 
and two-layer film are presented in A-B.l and A-B.2, respectively. In each section, 
the first program is to calculate a steady state solution as an input data of stability 
analysis. Unlike the program for steady state solution presented in Appendix A, the 
integration is done from the freeze line to the die. The second program is to calculate 
the eigenvalues of the system for stability analysis. The original length scale which is 
based on the radius of die, is changed to the one based on the freeze line height. The 
subroutines DLINRG and DGVLRG are from the IMSL.33 The program is used for 
both a boundary condition of fixed A and v^ (case 3) and a boundary condition of 
fixed B and v^- (case 1). If the absolute pressure in the program has a value other than 
zero, the program is for a case of fixed A and Vf boundary condition (case 3). If the 
absolute pressure is zero, the program is reduced to fixed B and Vf case (case 1). The 
program for case 4 (fixed A and T^) is easily obtained by letting T^’ zero and adding 
vj term instead. The subprogram FUNCTION RUNGE is not listed here, (see 
Appendix A) 

A-B. 1 Computer Program For Stability Analysis 
Of Newtonian Sinele-Laver Film 



C STEADY STATE SOLUTION OF NEWTONIAN FLUID SINGLE LAYER 
C BLOWN FILM EXTRUSION 
C BY FOURTH-ORDER RUNGE-KUTTA METHOD. 

C COMPUTATION FROM z=FL TO z=0 WITH DOUBLE PRECISION 
C THE SOLUTION IS USED AS INPUT DATA FOR STABILITY ANALYSIS 

C 



no 



Ill 



C DESCRIPTION OF VARIABLES ***• 

C 

C R(l)=r R(2)=r' R(3)=w 

C F(I)=r' F(2)=r" F(3)=w' 

C R(4)=aii R(5)=C33 
C 

IMPLICIT REAL*8 (A-H,0-Z) 

INTEGER RUNGE 
DIMENSION F(3),R(3) 

DATA DZ,ICOUNT,IFREQ/-5.D-4, 0,200/ 
OPEN(7,FILE='RNSS.DAT',STATUS=’UNKNOWN') 

C 

C **** PARAMETERS **** 

C 

C T= T-BAR, B= INFLATING PRESSURE, 

C FL=FREEZE LINE HEIGHT, BUR=BLOE UP RATIO, 

C TR= THICKNESS REDUCTION 
C 

PRINT *,'INPUT B,T,BUR,TR’ 

READ(*,*) B,T,BUR,TR 

TZ=T+B*BUR**2 

FL=5.D0 

WRITE (*,301) T,B,FL 
WRITE (7,*) T,B,FL 
C 

C INITIAL VALUES 
C 

3 Z=FL 
R(I)=BUR 
R(2)=0.D0 
R(3)=1.D0 

F(3)=-R(3)*((I.DO+R(2)**2)*(T+B*R(I)**2)/4.DO+R(2)/(R(I)*2.DO)) 

F(2)=(6.D0*R(2)+R(1)*(I.+R(2)*»2)*(T-3.D0*B*R(I)**2))/ 

&(2.D0*R(I)**2*(T+B*R(I)**2)) 

C 

1 CONTINUE 

WRITE(7,*) R(I),F(I),R(3),F(3),F(2) 

WRITE(*,102) Z,R(I),F(I),R(3),F(3),F(2) 

C 

C *♦♦♦ FOURTH-ORDER RUNGE-KLTTA METHOD 
C 

2 K=RUNGE(3,R,F,Z,DZ) 

IF(K.NE.l)GOTO 10 
F(1)=R(2) 

F(2)=(6.D0*R(2)+R(I)*(I.+R(2)**2)*(T-3.D0*B*R(I)**2))/ 

&(2.D0*R(1)**2*(T+B*R(1)**2)) 

F{3)=-R(3)*((1.D0+R(2)**2)*(T+B*R(I)**2)/4.D0+R(2)/(R(1)*2.D0)) 
GO TO 2 
C 

C ♦ ♦ * * PRINT INTERVAL * ♦ ♦ * 

C 

lOIF(Z.LT.-O.OI) GO TO 4 



112 



ICOUNT= ICOUNT+1 
IF (ICOUNT.NE.IFREQ) GO TO 2 
ICOUNT=0 
GOTO 1 
C 

C ***• TERMINATION ♦*** 

C 

4 CONTINUE 

WRITE(*,300) TZ,BUR,R(3),F(1) 

WRITE(7,*) TZ,BUR,R(3),F(1) 

STOP 

301 FORMAT(/3HT=,F6.3,4H B=,F6.3,5H FL=F6.3 
&//6X,lHX,9X,lHR,8X,7HR-PRIME,9X,lHW,8X,2HWP,8X,3HRDP) 
102 FORMAT (1H,F7.3,5D13.5) 

300 FORMAT(/4H TZ= F9.3,5H BUR= F6.3/ 

&4H TR=, F9.3,6H RPI=,F6.3) 

END 



The steady state solution computed above is given as input data in the 
following program to calculate eigenvalues. Terms from ALl to CRT represent the 
coefficients of the variables that appear in the linearized equations. The matrix 
coefficients A, B, C, AA, BB and CC are from the finite difference equations. 



C STABILITY ANALYSIS OF NEWTONIAN FLUID SINGLE-LAYER 
C BLOWN FILM EXTRUSION WITH DOUBLE PRECISION. 

C IF EXTERNAL PRESSURE P IS NOT ZERO, THIS SOLVES THE 
C CASE OF A FIXED A AND Vf . 

C IF P IS ZERO, THIS IS REDUCED TO THE CASE OF A FIXED B 
C ANDVf. 

C 

P ARAMETER(N=5 1 ,M=2 *N-2,M I =N- 1 , 
&LDA=M,LDB=M,LDCC=M1,LDCCI=MI) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/WORKSP/RWKSP 
DIMENSION RWKSP(178142) 

COMPLEX ALPHA(4*N-4) 

DIMENSION EVAL(4*N-4),BETA(2*N-2) 

DIMENSION R(N),W(N),RP(N),RDP(N),WP(N),SQ(N) 

DIMENSION A(LD A,M),B(LDB,M),C(M,M1), AA(M 1 ,M),BB(M 1 ,M), 
&CC(LDCC,M1),CCI(LDCCI,M1),CCIAA(MI,M),CCIBB(M1,M) 
DIMENSION BL1(N),BL2(N),BR1(N),BR2(N),BR3(N),BR4(N),BR5(N), 
&BR6(N),BR(7),CL1(N),CL2(N),CL3(N),CR1(N),CR2(N), 
&CR3(N),CR4(N),CR5(N),CR6(N),CR7(N),CR(8),AL1(N),AL2(N),AL3(N), 
&AR1(N),AR2(N),AR3(N),AR4(N),AR5(N),AR6(N) 

DIMENSION CCCIAA(M,M),CCCIBB(M,M) 



113 



DIMENSION RJM(M 1 ,M 1 ),CCIK(M 1 .M 1 ). AI(M 1 ,M 1 ), 
&RES(M1,M1),CCINEW(M1.M1).CCIN(M1.M1) 

EXTERNAL AMACH,DGVLRG,DLINRG 
CALL 1WKIN(178142) 

C 

C RNSS.DAT IS FROM THE ABOVE PRRGRAM FOR STEADY 
C STATE SOLUTION 
C 

OPEN(5,FlLE='RNSS.DAT',STATUS='OLD’) 

OPEN(9,FILE=’R.DAT',STATUS=’UNKNOWN') 

C 

C ********♦*♦**••**♦*♦♦♦♦♦♦•*♦**♦**♦♦♦♦♦♦♦♦* R F A n DATA 

C 

READ(5,*) T.BP.FL 
WRJTE(*,100) T.BP.FL 

lOOFORMATC T-BAR=’,F8.4,' BP=',F8.4,’ FL=',F8.4) 

DO I=N,1,-1 

READ(5,*) R(I),RP(I),W(I),WP(I),RDP(I) 

WRITE(*,1 1 1) I,R(I),RP(I),W(I),WP(I),RDP(I) 

111 FORMAT(I4,5F13.4) 

END DO 

READ(5,*) TZ,BUR,TR,RPI 
WRITE(*,*) TZ,BUR,TR,RPI 
H=l.DO/(DBLEfN-D) 

SF=1.D0/FL 

BP=BP/(SF**2) 

P=100.D0 
DO 2 J=1,N 

SQ(J)=DSQRT(1.D0+RP(J)**2) 

2 CONTINUE 
C 

C ******** CON\ERSION OF RNSS.DAT TO NEW SCALE 
C (SCALE FACTOR SF=1/FL) 

C 

DO 3 J=1,N 
R(J)=SF*R(J) 

W(J)=W(J)/TR 

WP(J)=WP(J)/(SF*TR) 

RDP(J)=RDP(J)/SF 

3 CONTINUE 
C 

Q **«********«:»« GIVING COEFFICIENTS OF THE MATRIX [A] TO [CJ 

DO 21 I=l,2*N-2 
DO 20 J=l,2*N-2 
A(I,J)=0.D0 

20 B(I,J)=0.D0 

21 CONTINUE 

DO 31 I=l,2*N-2 
DO 30 J=1,N-1 

30 C(U)=0.D0 

31 CONTINUE 
DO 41 1=1, N-1 
DO 40 J=l,2*N-2 



o n o 



114 



AA(1,J)=0.D0 

40 BB(U)=0.D0 

4 1 CONTINUE 
DO 51 I=1.N-1 
DO 50 J=1,N-1 
CC(I,J)=0.D0 

50 CCI(I,J)=0.D0 

51 CONTINUE 

**♦* GIVING VALUES OF AL-CR **** 

DO 55 J=1,N 
AL1(J)=-W(J)*SQ(J) 

AL2(J)=-R( J)* W(J) *RP( J)/SQ(J) 

AL3(J)=-R(J)*SQ(J) 

AR1(J)=SF*(-RP(J)/(R(J)*R(J))) 

AR2(J)=SF/R(J) 

AR3(J)=-SF*WP(J)/(W(J)*W(J)) 

AR4(J)=SF/W(J) 

AR5 ( J)=RP( J) * W( J)+R( J) * WP( J) 

AR6(J)=R(J)*W(J) 

C 

BL1(J)=2.D0*R(J)*RP(J)*W(J)/(SQ(J)**3) 

BL2(J)=-2.D0*R(J)/SQ(J) 

BR1(J)=(BP*R(J)*R(J)-T)/R(J) 

BR2(J)=2.D0*RP(J)*(T+BP*R(J)*R(J))/(SQ(J)*SQ(J)) 

BR3(J)=SF*(2.DOAV(J))*(WP(J)AV(J)+RP(J)/R(J))/(SQ(J)*SQ(J)) 

BR4(J)=SF*(2.D0AV(J))/(SQ(J)*SQ(J)) 

BR5(J)=2.D0*R(J)*WP(J)/(SQ(J)*SQ(J)) 

BR6(J)=-2.D0*R(J)*W(J)/(SQ(J)*SQ(J)) 

BR7(J)=R(J)**2-R(N)**2 

C 

TB=T+BP*R(J)*R(J) 

PAR1=-(RP(J)*W(J)-R(J)*WP(J))/{RP(J)*W(J)+2.D0*R(J)*WP(J)) 

PAR2=PAR1*TB/R(J) 

PAR3=-(RP(J)*W(J)+2.D0*R(J)*WP(J))/((R(J)*W(J))**2) 

CL 1 (J)=- W( J) ♦ RP( J)* RDP( J)/SQ( J) * * 3 

CL2(J)=W(J)*SQ(WJ)**2 

CL3 (J)=RDP(J)/SQ(J)-SQ(J)/R( J) 

CR1(J)=SF*(2.D0*RP(J)/R(J)-WP(J)/W(J))/R(J)**3 

CR2(J)=2.D0*BP*RP(J)-2.D0*RP(J)*RDP(J)/SQ(J)**4*SF/R(J) 

& ♦ (RP( J)/R( J)+2 .DO * WP( J)/W( J))-SF/R( J) ♦ * 3 
CR3(J)=-SF/R(J)*(RP(J)/R(J)+2.D0*WP(J)/W(J))/SQ(J)**2 
CR4(J)=-SF/R( J)AV( J) * (RDP( J) ♦ (R( J) * WP( J)+RP( J) * W( J))/ 
&(R(J)*W(J)*SQ(J)**2)+RP(J)/R(J)**2) 
CR5(J)=SF/(R(J)*W(J))*(-RDP(J)/SQ(J)**2+1.D0/R(J)) 

CR6(J)=-RDP( J) * WP( J)/SQ( J) * ♦ 2 -W( J) *RP( J)/R( J) * * 2 + WP( J)/R( J) 
CR7(J)=RDP( J) ♦ W( J)/SQ( J) * ♦ 2 
CR8(J)=SQ(J)**2 
55 CONTINUE 
C 

C ♦** COMPONENTS OF A, B & C FROM THE CONTINUITY EQ'N*** 
C 



DO 60 J=2,N-3 



non 



115 



A(J.J-1)=-AL2(J+1) 

A(J.J)=2.D0*H*AL1(J+1) 

A(J.J+1)=AL2(J+1) 

A(J,N- 1+J)=2.D0*H* AL3(J+ 1 ) 

C 

B(J,M)=-AR2(J+1) 

B(J,J)=2.D0*H*AR1(J+1) 

B(J.J+1)=AR2(J+1) 

B(J,N-2+J)=-AR4(J+l) 

B(J,N-l+J)=2.DO*H*AR3(J+l) 

B(J,N+J)=AR4(J+1) 

C 

C(J,M)=-AR6(J+1) 

C(J,J)=2.D0*H*AR5(J+1) 

C(J,J+1)=AR6(J+1) 

60 CONTINUE 
C 

J=1 

A(J,J)=2.D0*H*AL1(J+I) 

A(J,J+1)=AL2(J+1) 

A(J,N- 1+J)=2.D0*H* AL3(J+ 1 ) 

C 

B(J,J)=2.D0*H*AR1(J+1) 

B(J,J+1)=AR2(J+1) 

B(J,N-l+J)=2.DO*H*AR3(J+l) 

B(J,N+J)=AR4(J+1) 

C 

C(J,J)=2.D0*H*AR5(J+1) 

C(J,J+1)=AR6(J+1) 

C 

J=N-2 

A(J,J-1)=-AL2(J+1) 

A(J,J)=2.D0*H*AL1(J+1) 

A(J,J+1)=AL2(J+I) 

A(J,N-l+J)=2.DO*H*AL3(J+l) 

C 

B(J,J-1)=-AR2(J+1) 

B(J,J)=2.D0*H*AR1(J+1) 

B(J,J+1)=AR2(J+1) 

B(J,N-2+J)=-AR4(J+l) 

B(J,N-l+J)=2.DO*H*AR3(J+l) 

B(J,N+J)=AR4(J+1) 

C 

C(J,M)=-AR6(J+1) 

C(J,J)=2.D0*H*AR5(J+1) 

♦♦♦ WHEN K=N-1, BACK WARD DIFFERENCING IS USED 

A(N-l,N-2)=-AL2(N) 

A(N- 1 ,N- 1 )=H* AL I (N)+AL2(N) 

A(N-I,2*N-2)=H*AL3(N) 

C 

B(N-I,N-2)=-AR2(N) 

B(N- 1 ,N- 1 )=H* ARl (N)+AR2(N) 



non 



116 



B(N-l,2*N-3)=-AR4(N) 

B(N-1,2*N-2)=H*AR3(N)+AR4(N) 

C 

C(N-l,N-2)=-AR6(N) 

•** COMPONENTS OF A. B & C FROM THE AXIAL BALANCE *•* 



DO 70 J=2,N-3 
A(N-1+J,J-1)=-BL1(J+1) 
A(N-1+J,J+1)=BLI(J+1) 
A(N-l+J,N-l+J)=2.DO*H*BL2(J+l) 

C 

A(N- 1 + J,N-2)=2.D0* BL 1 (N) 

A(N- 1 + J,N- 1 )=-2 . DO * BL 1 (N) 
A(N-l+J,2*N-2)=-2.D0*H*BL2(N) 

C 

B(N-1+J,J-1)=-BR2(J+1) 

B(N- 1 +J, J)=2 ,D0*H*BR 1 (J+ 1 ) 
B(N-1+J,J+1)=BR2(J+1) 

B(N- 1 + J,N-2+J)=-BR4( J+ 1 ) 

B(N- 1 + J,N- 1 + J)=2 .DO *H*BR3 ( J+ 1 ) 
B(N-1+J,N+J)=BR4(J+1) 

C 

B(N-1+J,N-2)=2.D0*BR2(N) 

B(N- 1 +LN-1 )=-2.D0*(H*BR 1 (N)+BR2(N)) 
&-4.D0*P*H/R(N)*BR7(J+1) 
B(N-1+J,2*N-3)=2.D0*BR4(N) 
B(N-1+J,2*N-2)=-2.D0*(H*BR3(N)+BR4(N)) 

c 

C(N-1+J,M)=-BR6(J+1) 

C(N-l+J,J)=2.DO*H*BR5(J+l) 

C(N-1+J,J+1)=BR6(J+1) 

C 

C(N-1+J,N-2)=2.D0*BR6(N) 

70 CONTINUE 
C 

J=1 

A(N-1+J,J+1)=BL1(J+1) 

A(N-1+J,N-1+J)=2.D0*H*BL2(J+1) 

C 

A(N-1+J,N-2)=2.D0*BLI(N) 

A(N- 1 +J,N- 1 )=-2.D0 *BL 1 (N) 
A(N-1+J,2*N-2)=-2.D0*H*BL2(N) 

C 

B(N-1+J,J)=2.D0*H*BR1(J+1) 

B(N-1+J,J+1)=BR2(J+1) 

B(N-1+J,N-1+J)=2.D0*H*BR3(J+1) 

B(N-1+J,N+J)=BR4(J+1) 

c 

B(N-l+J,N-2)=2.DO*BR2(N) 

B(N-1+J,N-1)=-2.D0*(H*BR1(N)+BR2(N)) 

&-4.D0*P*H/R(N)*BR7(J+l) 

B(N-1+J,2*N-3)=2.D0*BR4(N) 

B(N-1+J,2*N-2)=-2.D0*(H*BR3(N)+BR4(N)) 



o o o o o r> 



117 



C 

C(N-1+J.J)=2.D0*H*BR5(J+1) 

C(N-1+J.J+I)=BR6(J+1) 

C 

C(N- 1 + J,N-2)=2. D0*BR6(N) 

C 

J=N-3 

A(N- 1 +J.N-2)=BL l(N-2)+2.D0*BL 1 (N) 
B(N-l+J,N-2)=BR2(N-2)+2.D0*BR2(N) 
B(N-l+J,2*N-3)=BR4(N-2)+2.D0*BR4(N) 
C(N-l+J.N-2)=BR6(N-2)+2.D0*BR6(N) 

C 

J=N-2 

A(N-1+J,M)=-BL1(J+1) 

A(N- 1 +J,N- 1 +J)=2.D0*H*BL2(J+ 1 ) 

C 

A(N-1+J,N-2)=2.D0*BL1(N) 

A(N- 1 + J,N- 1 )=BL 1 (N- 1 )-2.D0 *BL 1 (N) 
A(N-l+J,2*N-2)=-2.D0*H*BL2(N) 

C 

B(N-1+J,J-1)=-BR2(J+1) 

B(N- 1 +J.N-2+ J)=-BR4(J+ 1) 

C 

B(N-1 +J,N-2)=2.D0*(H*BR1 (N- 1 )+BR2(N)) 
B(N-1+J,N-1)=BR2(N-1)-2.D0*(H*BR1(N)+BR2(N)) 
&-4.D0*P*H/R(N)*BR7(J+l) 
B(N-1+J,2*N-3)=2.D0*(H*BR3(N-1)+BR4(N)) 
B(N-1+J,2*N-2)=BR4(N-1)-2.D0*(H*BR3(N)+BR4(N)) 

C 

C(N-1+J,J-1)=-BR6(J+1) 
C(N-1+J,N-2)=2.D0*(H*BR5(N-1)+BR6(N)) 

(2N-2)TH EQUATION FROM THE BOUNDARY CONDIION ♦** 



A(2*N-2,N-1)=H 

B(2*N-2,N-2)=l.DO/(R(N)*W(N)) 

B(2 ♦ N-2, N- 1 )=- 1 . D0/(R(N) ♦ W(N)) 

♦** COMPONENTS OF A A, BB & CC FROM THE RADIAL BALANCE 



DO 80 J=2,N-3 
AA(J,J-I)=-H*CLI(J+1) 
AA(J,J)=2.D0*H*H*CL2(J+I) 
AA(J,J+I)=H*CL1(J+1) 

AA(J,N- 1 +J)=2.D0*H*H*CL3(J+ 1 ) 

C 

BB(J, J- 1 )=2.D0*CR3 (J+ 1 )-H*CR2(J+ 1 ) 

BB(J,J)=2.D0*H*H*CRI(J+I)-4.D0*CR3(J+1) 

BB(J, J+ 1 )=H*CR2(J+ 1 )+2.D0*CR3(J+ 1 ) 

BB(J,N- 1 )=-4.D0*H* *2*P/R(N)*CR8( J+ 1 ) 

BB(J,N-2+J)=-H*CR5(J+I) 

BB(J,N-I+J)=2.D0*H*H*CR4(J+I) 

BB(J,N+J)=H*CR5(J+I) 



oooo ooo oon 



118 



CC(J,M)=-H*CR7(J+1) 

CC(J,J)=2.D0*H*H*CR6(J+1) 

CC(J,J+l)=H*CR7(J+l) 

80 CONTINUE 
C 

J=1 

AA(J,J)=2.D0*H*H*CL2(J+1) 

AA(J,J+l)=H*CLl(J+l) 

AA(J,N-l+J)=2.DO*H*H*CL3(J+l) 

C 

BB(J,J)=2.DO*H*H*CR1(J+1)-4.DO*CR3(J+1) 

BB(J,J+1)=H*CR2(J+1)+2.D0*CR3(J+1) 

BB(J,N-1)=-4.D0*H**2*P/R(N)*CR8(J+1) 

BB(J,N-1+J)=2.D0*H*H*CR4(J+1) 

BB(J.N+J)=H*CR5(J+1) 

C 

CC(J,J)=2.D0*H*H*CR6(J+1) 

CC(J,J+1)=H*CR7(J+1) 

C 

J=N-2 

AA(J,J-1)=-H*CL1(J+1) 

AA(J,J)=2.D0*H*H*CL2(J+1) 

AA(JJ+1)=H*CL1(J+1) 

AA(J,N-l+J)=2.DO*H*H*CL3(J+l) 

C 

BB(J,J-1)=2.D0*CR3(J+1)-H*CR2(J+1) 

BB(J,J)=2.D0*H*H*CR1(J+1)-4.D0*CR3(J+1) 

BB(J,J+1)=H*CR2(J+1)+2.D0*CR3(J+1) 

&-4.D0*H**2*P/R(N)*CR8(J+l) 

BB(J,N-2+J)=-H*CR5(J+l) 

BB(J,N-1+J)=2.D0*H*H*CR4(J+1) 

BB(J,N+J)=H*CR5(J+1) 

C 

CC(J,J-1)=-H*CR7(J+1) 
CC(J,J)=2.D0*H*H*CR6(J+I) 

*•* (N-1) TH EQUATION FROM dTz EXPRESSION 

AA(N-I,2*N-2)=2.*H*R(N) 

BB(N- 1 ,N- 1 )=-4. ♦ SF*H* WP(N)/(R(N) * W(N)) 
BB(N-I,2*N-3)=2.*SF/W(N) 
BB(N-I,2*N-2)=-(2.*SF/W(N))*(H*WP(N)/W(N)-l.) 
CC(N-I,N-2)=-2.*R(N)*W(N) 

CC(N-I,N-I)=-H 






INVERSE OF [CC] 



CALL DLINRG(MI,CC,LDCC,CCLLDCCI) 

RECURSION FORMULA WHEN THE INVERSE IS VERY ILL- 
CONDITIONED 

ITER=0 

DO 300 I=I,N-I 



119 



DO 300 J=l,N-l 
RIM(I,J)=0. 

300 IF(I.EQ.J) R1M(I,J)=1. 

GO TO 499 
299 ITER=ITER+1 
WRITE(*,*) ITER 
C 

DO 400 1=1, N-1 
DO 400 J=1,N-1 
CCIK(I,J)=0. 

DO410K=l,N-l 

410 CCIK(I, J)=CCIK(I,J)+CC(I,K)*CCI(K. J) 

400 CONTINUE 
C 

DO 500 1=1, N-1 
DO 500 J=1,N-1 

500 AI(I,J)=2.*RIM(I,J)-CCIK(I,J) 

C 

DO 600 1=1, N-1 
DO 600 J=1,N-1 
CCIN(I,J)=0. 

DO610K=l,N-l 

610 CCIN(I, J)=CCIN(I, J)+CCI(I,K)* AI(K, J) 

600 CONTINUE 

C 

DO 611 I=1,N-1 
D0 612 J=1,N-1 
612 CCI(I,J)=CCIN(I,J) 

6 1 1 CONTINUE 
C 

499 DO 622 I=1,N-1 
DO 622 J=1,N-1 
CCINEW(I,J)=0. 

DO 623 K=1,N-1 

623 CCINEW(I,J)=CCINEW(I,J)+CC(I,K)*CCI(K,J) 

622 CONTINUE 
C 

DO 650 1=1, N-1 
DO 650 J=1,N-1 

RES(I,J)=RIM(I,J)-CCINEW(I,J) 

IF (DABS(RES(I,J)).GT.l.D-8) GO TO 299 
650 CONTINUE 

WRITE(*,911) ITER 

911 FORMAT(lX,'NO. OF ITERATION FOR INVERSE= ’,13) 
C 

Q ♦♦*♦♦♦*♦»♦♦♦♦♦♦♦♦******* MULTIPLYING [CCI][AA] 
C & [CCI][BB] 

DO 700 1=1, N-1 
DO 700 J=l,2*N-2 
CCIAA(I,J)=0. 

CCIBB(I,J)=0. 

DO710K=l,N-l 

CCIAA(I,J)=CCIAA(I,J)+CCI(I,K)*AA(K,J) 

710CCIBB(I,J)=CCIBB(I,J)+CCI(I,K)*BB(K,J) 



120 



700 CONTINUE 
C 

C 



COMPUTE [CIlCCIAAj AND (C|(CCIBB] 



DO 830 I=l,2*N-2 
DO 830 J=l,2*N-2 
CCCIAA(I,J)=0. 

CCCIBB(I,J)=0. 

DO 831 K=1,N-1 

CCCIAA(U)=CCCIAA(I,J)+C(I,K)‘CCIAA(K,J) 

83 1 CCCIBB(I,J)=CCCIBB(I,J)+C(I,K)*CCIBB(K,J) 

830 CONTINUE 
C 

Q ♦**♦♦*♦**♦♦**♦*♦♦*♦♦** 1A]=[A]-[CCCIAA] AND {B]=[BJ-rCCCIBBl 
C 



DO 860 I=l,2*N-2 
DO 861 J=l,2*N-2 
A(I,J)=A(I,J)-CCCIAA(I,J) 
861 B(U)=B(I,J)-CCCIBB(U) 



860 CONTINUE 
C 

C 



COMPUTE EIGENVALUE FROM DG\TRG 



CALL DGVLRG(M,B,LDB,A,LDA.ALPHA,BETA) 

DO 71 I=2*N-2,1,-1 
IF(BETA(I).NE.0.0) THEN 
EVAL(2*I-1)=DBLE(ALPHA(2*I-1))/BETA(I) 
EVAL(2*I)=DBLE(ALPHA(2*I))/BETA(I) 

ELSE 

E VAL(2*I- 1 )= AMACH(2) 

EVAL(2*I)=AMACH(2) 

ENDIF 

71 CONTINUE 
C 

Q *♦♦*♦**♦♦♦♦♦*♦*♦♦♦«♦♦**♦*** print THE EIGENVALUES 
C 



WRITE(*,901) N,2*N-2 
WRITE(9,90I) N,2*N-2 

901 FORMAT(//5X,'NUMBER OF NODE= ’,I3//5X,'MATRIX SIZE= ’,13) 
BP=BP*SF*SF 

WRITE(*,902) TZ,BP,BUR,TR.T.RPI 
WRITE(9,902) TZ,BP,BUR,TR.T,RPI 

902 FORMAT(//lX,'DRAW FORCE= ',F9.3,5X,’B= ',F6.3/1X,'BUR= 
&F6.3,15X,TR= ',F9.3/1X,T= ',F8.4,15X,'RPI= ’,F10.6,/) 

WRITE(*,985) P,FL 
WRITE(9,985) P,FL 

985 FORMAT(2X,’PRESSURE= ',F6. l,2x.'FL= ',f4.1) 

WRITE(*,986) 

WRITE(9,986) 

986 FORMAT(//,3X,’NO.',13X,'RE.AL',16X,’IMAGV) 

K=0 

DO J=2*N-2,2*N-7,-l 
K=K+1 

EVR=EVAL(2*J-1) 



121 



EVI=EVAL(2*J) 
WRITE(*,987) K.EVR.EVI 
END DO 

987 FORMAT(3X,I2.2X,2E20.9) 
STOP 
END 



A-B.2 Computer Program For Stability Anavsis 

Of Two-Laver Film 



C STEADY SOLUTION OF TWO-LAYER FILM 
C BY FOURTH-ORDER RUNGE -KUTTA METHOD 
C AS INPUT DATA FOR STABILITY ANALYSIS 
C USING NEW SCALE FACTOR SF=Ro/L 
C INTEGRATION IS FROM z=0 TO z=FL 
^ ♦♦♦♦*♦**♦♦♦♦♦♦♦♦*♦♦♦♦♦♦♦♦♦♦♦*♦♦**♦♦♦♦*♦*♦*♦************** 
C 

IMPLICIT REAL*8 (A-H,0-Z) 

INTEGER RUNGE 
DIMENSION F(6),R(6) 

ICOUNT=0 

DZ=.0001 

IFREQ=250 

OPEN(4,FILE='COEX.DAr,STATUS='UNKNOWN’) 

C 

C ♦*** DESCRIPTION OF VARIABLES **♦* 

C 

C R(l)=r R(2)=r' R(3)=w 

C R(4)=SIGI1 R(5)=SIG33 R(6)=TAU22 

C F(I)=r' F(2)=r" F(3)=w' 

C F(4)=SIGir F(5)=SIG33’ F(6)=TAU22‘ 

C 

C 

C PARAMETERS 

C 

C T= T-BAR, B= INFLATING PRESSURE, 

C D=RELAXATION TIME, FR=FLOW RATE RATIO, 

C FL=FREEZE LINE HEIGHT 
C 

C INITIAL GUESS FOR R-PRIME ♦*** 

C 

READ (♦,*)D,FR,B,T,RPI 

VISCO=I.DO 

FL=5.D0 

G=( 1 . DO-FR) * VIS CO/FR 
SF=I./FL 

WRITE (4,*) T,B,D,FR,VISCO,FL 
IP=0 



o o 



122 



INITIAL VALUES 
3 Z=O.DO 
R(1)=SF 
R(2)=RPI 
R(3)=1.D0 
R(4)=-. I 
R(5)=.l 
R(6)=0. 

Tl=R(l)/SF*R(3)*DSQRT(l.DO+R(2)*R(2))/D 
T2=R(2)/R(1)*SF 

T3I=-T2/2.D0+T1*D*R(4)/(2.D0*G) 
T3=T3I-(1.D0+R(2)*R(2))*(T+B*(R(1)/SF)**2)/(4.D0*FR*G) 
F(1)=R(2) 

F2I=(R(5)+G*(T2-T3)/(T1*D))/(R(4)-G*(T2+2.D0*T3)/(T1*D)) 
F(2)=(1.DO+R(2)*R(2))*(-2.DO*B*R(1)/SF 
&/(T+B ♦ (R( 1 )/SF) ♦ * 2 )+F2 1/R( 1 ) * SF)/S F 
F(3)=T3*R(3)/SF 

F4I=-(T1+2.D0*(T2+T3))*R(4)-2.D0*(2.D0*T3+T2)*R(6) 
F(4)=(F4I-(2.D0*T3+T2)/D)/SF 

F(5)=((2.D0*T2-T1)*R(5)+(T2-T3)*(2.D0*R(6)+I.D0/D))/SF 
F(6)=(T3/D+(2.D0*T3-T1)*R(6))/SF 

1 CONTINUE 

WRITE(4,*) R(I),R(3),R(4),R(5),R(6), 

& F(1),F(2),F(3),F(4),F(5),F(6) 

IF(Z.GE.1.D0-.1D0*DZ) GO TO 4 
C 

C FOURTH ORDER RUNGE-KUTTA 

2 K=RUNGE(6,R,F,Z,DZ) 

IF (K.NE.l) GO TO 10 

Tl=R(l)/SF*R(3)*DSQRT(l.DO+R(2)*R(2))/D 

T2=R(2)/R(1)*SF 

T3I=-T2/2.D0+T1*D*R(4)/(2.D0*G) 

T3=T3I-(1.D0+R(2)*R(2))*(T+B*(R(1)/SF)**2)/(4.D0*FR*G) 

F(1)=R(2) 

F2I=(R(5)+G*(T2-T3)/(T1*D))/(R(4)-G*(T2+2.D0*T3)/(T1*D)) 
F(2)=( 1 .D0+R(2)*R(2))*(-2.D0*B*R( 1 )/SF 
&/(T+B *(R( 1 )/SF)* * 2)+F2I/R( 1 ) ♦ SF)/SF 
F(3)=T3*R(3)/SF 

F4I=-(T1+2.D0*(T2+T3))*R(4)-2.D0*(2.D0*T3+T2)*R(6) 

F(4)=(F4I-(2.D0*T3+T2)/D)/SF 

F(5)=((2.D0*T2-T1)*R(5)+(T2-T3)*(2.D0*R(6)+1.D0/D))/SF 
F(6)=(T3/D+(2.D0*T3-T1)*R(6))/SF 
GO TO 2 
C 

C PRINT INTERVAL ♦*** 

10 CONTINUE 

ICOUNT= ICOUNT+1 
IF (ICOUNT.NE.IFREQ) GO TO 2 
ICOUNT=0 
GOTO 1 
C 

C **♦* TERMINATION 



123 



4 IF(IP.EQ.O) GO TO 5 
RPF=R(2) 

IF (DABS(RPF).LE.l.D-2) GO TO 5 
WRITE (*,103) RPI.RPF 
RPI=RPI+l.D-2 
GO TO 3 

5 CONTINUE 
TZ=T+B*(R(1)/SF)**2 
TR=I.D0/R(3) 

BUR=R(1)*FL 

WRITE (4,*) TZ,BUR,TR,RPI 
WRITE(*,221) TZ,B,BUR,TR 

221 FORMAT(/13H DRAW FORCE= ,F9.3,5X,4H B= ,F6.3/ 

&6H BUR= ,F6.3, 12X,5H Tr= ,F9.3/) 

STOP 

101 FORMAT (/5H DZ= ,F7.5/2X,3H X ,8X,1HR,7X.7HR-PRIME, 
&6X,1HW,8X,5HSIG11,6X,5HSIG33,6X.5HTAU22) 

102 FORMAT (I7.3,5(lx,dl0.3),lx,d9.3) 

103 FORMAT (6H RPI=,D13.5,6H RPF=.D13.5) 

104 FORMAT(/7H T-BAR=, F9.3,5H B=.F6.3,9H LAMDA=,F6.3/ 

6 9HF-RATIO=,F6.3,10H V-RATIO=. F6.3, 13H FROST-LINE= F6.3, ) 
END 



The above solution is used as an input data in the following program. Since the 

parameters are rescaled in the above program, there is no modification in the following 
program. 



C STABILITY ANALYSIS OF THE TWO-LAYER 
C BLOWN FILM EXTRUSION AT A FIXED A AND Vf 
C DOUBLE PRECISION 

C ***************’*‘**********^^**** *♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦* 

c 

PARAMETER(N=41) 

COMMON/WORKSP/RWKSP 
REAL RWKSP(361522) 

COMPLEX ALPHA(10*N-8) 

REAL*8 EVAL(10*N-8),BETA(5*N-4) 

REAL*8 AMACH 

EXTERNAL AMACH,DGVLRG,DLINRG 
REAL*8 R(N),W(N),RP(N),RDP(N),WT(N), 
&S11(N),S33(N),TAU(N),S11P(N),S33P(N),TAUP(N), 
&SQ(N),A(5*N-4,5*N-4),B(5*N-4,5*N-4),C(5*N-4,N-1), 
&AA(N-1,5*N-4),BB(N-1,5*N-4),CC(N-1,N-1),BL1(N),BL2(N), 
&BR1(N),BR2(N),BR3(N),BR4(N),BR5(N),BR6(N),BR7(N),BR8(N) 

REAL*8 CR10(N),CL1(N),CL2(N),CL3(N),CR1(N),CR2(N),CR3(N).CR4(N), 
&CR5(N),CR6(N),CR7(N),CR8(N),CR9(N),CR10(N),AL1(N),AL2(N), 
&AL3(N),AR1(N),AR2(N),AR3(N),AR4(N),AR5(N),AR6(N), 



124 



«feDLl(N),DL2(N),DL3(N).DRl(N).DR2(N),DR3{N).DR4(N),DR5(N) 
REAL*8 DR6(N),DR7(N),DR8(N).EL1(N),EL2(N),EL3(N), 
&ER1(N),ER2(N),ER3(N),ER4(N),ER5(N).ER6(N).ER7(N),ER8(N). 
&FL1(N),FL2(N),FR1(N),FR2(N),FR3(N),FR4(N),FR5(N).FR6(N), 

&CCI AA(N- 1 ,5 *N-4),CCCIBB(5 *N-4,5 •N-4).CCI(N- 1 .N- 1 ) 

REAL*8 CCIBB(N-l,5*N-4),CCCIAA(5*N-4,5*N-4), 
&RIM(N-1,N-1),CCIK(N-1,N-1),AI(N-1,N-1), 

&RES(N- 1 ,N- 1 ),CCINE W(N- 1 ,N- 1 ), CCIN(N- 1 , N- 1 ) 

REAL*8 T,BP,D,F.VISC0,FL,TZ,BUR,TR,RPI,H,SF.RF,P,TB.VV 
CALL IWKIN(361522) 

OPEN(5,FILE='C2.DAT',STATUS=’OLD’) 

0PEN(9,FILE='TLS1.DAT',STATUS='UNKN0WN') 

C 

C *’*******’*’*********♦♦♦♦♦♦♦♦*♦♦♦♦♦♦*♦♦♦♦♦♦♦* READ INPUT DATA 
C 



READ(5,*) T,BP,D,F,VISCO,FL 
DO 1=1, N 

READ(5,») R(I), W(I),S 1 1 (I),S33(I),TAU(I), 

& RP(I),RDP(I),WP(I),S1 1P(I),S33P(I).TAUP(I) 

WR1TE(*,*) I,R(I),RP(I),W(I) 

END DO 

READ(5,*) TZ,BUR,TR,RPI 
H=l.DO/(DBLE(N-l)) 

SF=1.D0/FL 

BP=BP/(SF**2) 

P=100.D0 
DO 2 J=1,N 

2 SQ(J)=DSQRT(1.D0+RP(J)**2) 

C 

Q ♦**♦♦******♦,* Qjyg COEFFICIENTS OF THE MATRIX [A] TO [C] 

DO 21 I=l,5*N-4 
DO 20 J=l,5*N-4 
A(I,J)=0. 

20 B(I,J)=0. 

21 CONTINUE 
DO 31 J=1,N-1 
DO 30 I=l,5*N-4 

30 C(I,J)=0. 

31 CONTINUE 

DO 40 l=l,5*N-4 
DO 41 J=1,N-1 
BB(J,I)=0. 

41 AA(J,I)=0. 

40 CONTINUE 
DO 51 1=1,N-1 
DO 50 J=1,N-1 
CC1(1,J)=0. 

50 CC(1,J)=0. 

51 CONTINUE 
C 

C GIVING VALUES OF AL-FR ♦♦♦♦*♦*♦♦♦*♦♦♦♦**♦*♦♦*♦ 

C 



125 



DO 55 J=1,N 
AL1(J)=-W(J)*SQ(J) 

AL2(J)=-R( J) * W( J) *RP( J)/SQ( J) 

AL3(J)=-R(J)*SQ(J) 

AR 1 ( J)=-RP( J)/R( J) ♦ * 2 • SF 
AR2(J)=SF/R(J) 

AR3 ( J)=-SF * WP( J)/W( J) **2 
AR4(J)=SF/W(J) 

AR5(J)=RP(J)*W(J)+R(J)*WP(J) 

AR6(J)=R(J)*W(J) 

BL1(J)=(1.-F)*2.D0*R(J)*RP(J)*W(J)/SQ(J)**3 
BL2(J)=-( 1 . -F)*2.D0*R(J)/SQ(J) 
BR1(J)=-2./SF*F*W(J)*S11(J)/SQ(J)+2.*SF*(1.-F)/SQ(J)**2 
&/R(J)*(RP(J)/R(J)+2.*WP(J)/W(J))+2.*BP*R(J) 
BR2(J)=2./SF*F*R(J)*W(J)*S11(J)*RP(J)/SQ(J)**3 
&-4.*SF*(l.-F)*RP(J)/SQ(J)**4*(RP(J)/R(J)+2.*WP(J)/W(J)) 
BR3(J)=-2./SF*F*Sll(J)*R(J)/SQ(J)+2.*SF*(l.-F)* 
&(R(J)*WP(J)+RP(J)*W(J))/(R(J)*W(J)**2*SQ(J)**2) 
BR4(J)=2.*SF*(1.-F)/SQ(J)**2/W(J) 

BR5(J)=2 . * ( 1 . -F)*R( J)* WP( J)/SQ(J) * * 2 
BR6(J)=-2.*(1.-F)*R(J)*W(J)/SQ(J)**2 
BR7(J)=-2./SF*F*R(J)*W(J)/SQ(J) 

BR8(J)=R(J)**2 

1 

CL 1 ( J)=( 1 . -F) * SQ( J) ♦ W( J)/R( J) ♦ *2 
CL2(J)=-( 1 . -F)*RP(J)*RDP(J)* W(J)/SQ(J)* *3 
CL3(J)=(1.-F)*(-SQ(J)/R(J)+RDP(J)/SQ(J)) 
CR1(J)=F/SF*W(J)*S33(J)*SQ(J)/R(J)**2 
&+SF*( 1 . -F) * (2. ♦RP(J)/R(J)- WP( J)/W( J))/R( J) * * 3 
CR2(J)=2.*BP*RP(J)-F/SF*W(J)*RP(J)*(RDP(J)*S11(J)/SQ(J)**3 
&+S33(J)/R(J))+SF*(l.-F)*2.*RP(J)*RDP(J)/R(J)/SQ(J)**4 
& *(RP(J)/R(J)+2. ♦ WP(J)/W( J))-SF*( 1 . -F)/R( J) ♦ * 3 

CR3(J)=F/SF*W(J)*S11(J)/SQ(J)-SF*(1.-F)/SQ(J)**2*(RP(J)/'R(J)+ 

&2.*WP(J)/W(J))/R(J) 

CR4(J)=F/SF*(RDP(J)*S11(J)/SQ(J)-S33(J)*SQ(J)/R(J)) 

&-SF*(l.-F)/(R(J)*W(J))*(RDP(J)*(RP(J)*W(J)+R(J)*WP(J))/ 

&(R(J)*W(J)*SQ(J)**2)+RP(J)/R(J)**2) 

CR5(J)=SF*(1.-F)*(-RDP(J)/(R(J)*W(J)*SQ(J)**2) 

&+l./(R(J)**2*W(J))) 

CR6(J)=(1.-F)*(-RDP(J)*WP(J)/SQ(J)**2-W(J)*RP(J) 

&/R(J)**2+WP(J)/R(J)) 

CR7( J)=( 1 . -F) ♦ RDP( J) ♦ W(J)/SQ( J) * * 2 
CR8(J)=F/SF*W(J)^RDP(J)/SQ(J) 

CR9(J)=-F/SF*W(J)*SQ(J)/R(J) 

CR10(J)=SQ(J)**2 

DL 1 (J)=SF*RP(J)/SQ(J)* *2*(2. *D*S 1 1 (J)+2. *D*TAU(J)+ 1 .) 
DL2(J)=-SF* (2. *D *TAU(J)+ 1 . )/W( J) 

DL3(J)=-D*SF 

DR1(J)=-SF**2*D/(R(J)*W(J)*SQ(J)**3)*(RP(J)*S11P(J)+2.*(S11(J)+ 

&TAU(J))*RP(J)*(WP(J)/W(J)+RP(J)/R(J))+2.*TAU(J)*WP(J)*RP(J) 

&AV(J))-SF**2*RP(J)/SQ(J)**3/(R(J)*W(J))*(RP(J)/R(J)+2.D0* 

&WP(J)/W(J)) 



non 



126 



DR2(J)=-SF**2*WP(J)*(2.*D*TAU(J)+I.)/(R(J)*W(J)**3*SQ(J)) 

DR3(J)=SF**2*((2.*D*TAU(J)+1.)/(R(J)*W(J)**2*SQ(J))) 

DR4(J)=SF*D/SQ(J)*(S11P(J)+2.*TAU(J)*WP(J)/W(J))+WP(J) 

&/W(J)/SQ(J)*SF 

DR5(J)=-SF*(2. *D*S 1 1 (J)+2.*D*TAU(J)+ 1 ,)/SQ(J) 
DR6(J)=1.+2.*SF**2*D/SQ(J)*(R(J)*WP(J)+RP(J)*W(J)) 
&/(R(J)*W(J))**2 

DR7(J)=SF**2*D/(R(J)*W(J)*SQ(J)) 

DR8(J)=2 . * SF* *2 *D *(RP(J)/R( J)+2 . • WP(J)/W( J))/(R( J)* W( J)* SQ( J)) 
C 

EL1(J)=SF*((2.*D*S33(J)+2.*D*TAU(J)+1.)/R(J)) 

EL2(J)=-SF*(2.*D*TAU(J)+1.)AV(J) 

EL3(J)=-SF*D 

ER1(J)=SF**2*(RP(J)/(R(J)**3*W(J)*SQ(J))*(2.*D*S33(J)+2.*D* 

«&TAU(J)+1.)) 

ER2(J)=-SF**2*RP(J)/(R(J)*W(J)*SQ(J)**3)*(D*S33P(J)-2.*D*RP(J) 

&/R(J)*(S33(J)+TAU(J))+2.*D*TAU(J)*WP(J)/W(J)-RP(J)/R(J) 

&+WP(J)/W(J))-(2.*D*S33(J)+2.*D*TAU(J)+l.)/ 

&(R( J) * * 2 ♦ W( J) * SQ( J)) * SF* ♦ 2 

ER3(J)=-SF**2*WP(J)/(R(J)*W(J)**3*SQ(J))*(2.*D*TAU(J)+1.) 

ER4(J)=SF**2*((2.*D*TAU(J)+1.)/(R(J)*W(J)**2*SQ(J))) 

ER5(J)=SF*D/SQ(J)*(S33P(J)-2.*S33(J)*RP(J)/R(J)-2.*TAU(J)* 

&(RP(J)/R(J)-WP(J)/W(J)))-SF*(RP(J)/R(J)-WP(J)/W(J))/SQ(J) 

ER6( J)= 1 .-2 . ♦ SF* ♦ 2 *D *RP( J)/(R( J) * * 2 * W( J) * SQ( J)) 

ER7( J)=SF* ♦ 2 • D/(R( J) * W( J) ♦ SQ( J)) 

ER8(J)=-2.*SF**2*D*(RP(J)/R(J)-WP(J)/W(J))/(R(J)*W(J)*SQ(J)) 

C 

FL1(J)=SF*(2.*D*TAU(J)+1.)/W(J) 

FL2(J)=-SF*D 

FR1(J)=SF**2*RP(J)/(R(J)*W(J)*SQ(J)**3)*(-D*TAUP(J)+ 

&2.*D*TAU(J)*WP(J)/W(J)+WP(J)AV(J)) 

FR2(J)=SF*(D*TAUP(J)-WP(J)AV(J)*(2.*D*TAU(J)+1.))/SQ(J) 

FR3(J)=SF**2*WP(J)/(R(J)*W(J)**3*SQ(J))*(2.*D*TAU(J)+1.) 

FR4(J)=-SF**2*(2.*D*TAU(J)+1.)/(R(J)*W(J)**2*SQ(J)) 

FR5(J)=1.-2.*SF**2'*D*WP(J)/(R(J)*W(J)**2*SQ(J)) 

FR6( J)=SF* ♦ 2 * D/(R( J) * W( J) * SQ( J)) 

55 CONTINUE 

♦♦♦ COMPONENTS OF A,B,C FROM CONTINUITY EQ’N ♦♦♦ 



DO 60 J=2,N-3 
A(J,M)=-AL2(J+1) 
A(J,J)=2.D0*H*ALI(J+I) 
A(J,J+1)=AL2(J+I) 

A(J,N- 1 +J)=2. DO *H* AL3 (J+ 1 ) 
C 

B(J,J-I)=-AR2(J+I) 

B(J,J)=2.D0*H*ARI(J+I) 

B(J,J+1)=AR2(J+I) 

B(J,N-2+J)=-AR4(J+I) 

B( J,N- 1 + J)=2 .DO ♦ H * AR3 ( J+ 1 ) 
B(J,N+J)=AR4(J+1) 



C(J,M)=-AR6(J+1) 



nnooo o non 



127 



C(J,J)=2.DO*H*AR5(J+l) 

C(J,J+1)=AR6(J+1) 

60 CONTINUE 
C 

J=1 

A(J,J)=2.D0*H*ALI(J+1) 

A(J.J+1)=AL2(J+1) 

A(J,N-l+J)=2.DO*H*AL3(J+l) 

C 

B(J,J)=2.D0*H*AR1(J+1) 

B(J,J+1)=AR2(J+1) 

B(J,N-1+J)=2.D0*H*AR3(J+1) 

B(J,N+J)=AR4(J+1) 

C 

C(J,J)=2.DO*H*AR5(J+l) 

C(J,J+1)=AR6(J+1) 

C 

J=N-2 

A(J,J-1)=-AL2(J+1) 

A(J,J)=2.D0*H*AL1(J+1) 

A(J,J+1)=AL2(J+1) 

A(J,N-l+J)=2.DO*H*AL3(J+l) 

C 

B(J,M)=-AR2(J+1) 

B(J,J)=2.D0*H*AR1(J+I) 

B(J,J+1)=AR2(J+1) 

B(J,N-2+J)=-AR4(J+l) 

B(J,N-1+J)=2.D0*H*AR3(J+1) 

B(J,N+J)=AR4(J+1) 

C 

C(J,J-1)=-AR6(J+1) 

C(J,J)=2.DO*H*AR5(J+l) 

WHEN J=N-1, BACKWARD DIFFERENCE IS USED 

A(N-I,N-2)=-AL2(N) 

A(N- 1 ,N- 1 )=H* AL 1 (N)+ AL2(N) 
A(N-I,2*N-2)=H*AL3(N) 



B(N-I,N-2)=-AR2(N) 

B(N- 1 ,N- 1 )=H* AR 1 (N)+ AR2(N) 

B(N-I,2*N-3)=-AR4(N) 

B(N-I,2*N-2)=H*AR3(N)+AR4(N) 



C(N-l,N-2)=-AR6(N) 



*** COMPONENTS OF A,B & C FROM AXIAL FORCE BALANCE 

DO 61 J=2,N-4 
A(N-1+J,J+1)=BL1(J+1) 

A(N-1+J,J-1)=-BL1(J+1) 

A(N-1+J,N-2)=2.D0*BL1(N) 

A(N- 1 + J.N- 1 )=-2 . D0*BL 1 (N) 



o o 



128 



A(N-1+J,N-1+J)=2.D0*H*BL2(J+1) 

A(N-1+J,2*N-2)=-2.D0*H*BL2(N) 

C 

B(N-1+J,J+1)=BR2(J+1) 

B(N-1+J,J)=2.D0*H’*BR1(J+1) 

B(N-1+J,M)=-BR2(J+1) 

B(N- 1 +J,N-2)=2.D0*BR2(N) 

B(N- 1 +J,N- 1 )=-2.D0*H*(BR 1 (N)-2 . *P/R(N)*(BR8(N)-BR8(J+ 1 ) 
&))-2.D0*BR2(N) 

B(N-1+J,N+J)=BR4(J+1) 

B(N-1+J,N-1+J)=2.D0*H*BR3(J+1) 

B(N- 1 +J,N-2+J)=-BR4(J+ 1 ) 

B(N- 1 + J,2 *N-3 )=2 . D0*BR4(N) 
B(N-l+J,2*N-2)=-2.D0*H*BR3(N)-2.D0*BR4(N) 
B(N-1+J,2*N-1+J)=2.D0*H*BR7(J+1) 
B(N-1+J,3*N-2)=-2.D0*H*BR7(N) 

C 

C(N-1+J,J+1)=BR6(J+I) 

C(N-1+J,J)=2.D0*H*BR5(J+1) 

C(N-1+J,J-1)=-BR6(J+1) 

C(N-1+J,N-2)=2.D0*BR6(N) 

61 CONTINUE 



AT J=1 
J=1 

A(N-1+J,J+1)=BL1(J+1) 

A(N-1 +J,N-2)=2.DO*BL 1 (N) 

A(N- 1 + J,N- 1 )=-2. DO *BL 1 (N) 
A(N-1+J,N-1+J)=2.D0*H*BL2(J+1) 
A(N-1+J,2*N-2)=-2.D0*H*BL2(N) 

C 

B(N-1+J,J+1)=BR2(J+1) 

B(N-1+J,J)=2.D0*H*BR1(J+1) 

B(N-1+J,N-2)=2.D0*BR2(N) 

B(N-1+J,N-1)=-2.D0*H*(BR1(N)-2.*P/R(N)*(BR8(N)-BR8(J+1) 

&))-2.D0*BR2(N) 

B(N-1+J,N+J)=BR4(J+1) 

B(N-1+J,N-1+J)=2.D0*H*BR3(J+1) 

B(N-1+J,2*N-3)=2.D0*BR4(N) 

B(N-1+J,2*N-2)=-2.D0*H*BR3(N)-2.D0*BR4(N) 

B(N-1 +J,2*N-1+J)=2.D0*H*BR7(J+ 1) 
B(N-l+J,3*N-2)=-2.D0*H*BR7(N) 

C 

C(N-1+J,J+1)=BR6(J+1) 

C(N-1+J,J)=2.D0*H*BR5(J+1) 

C(N-l+J,N-2)=2.DO*BR6(N) 

C 

C AT J=N-3 
J=N-3 

A(N- 1+J, J+ 1)=BL 1 (J+ l)+2.DO*BL 1 (N) 

A(N-1+J,J-1)=-BL1(J+1) 

A(N-1 +J.N- 1 )=-2.D0*BL 1(N) 

A(N-1+J,N-1+J)=2.D0*H*BL2(J+1) 

A(N-l+J,2*N-2)=-2.D0*H*BL2(N) 



noonn non 



129 



C 

B(N-1 +J.J+ 1)=BR2(J+ 1 )+2.D0*BR2(N) 
B(N-1+J,J)=2.D0*H*BR1(J+1) 

B(N-1+J,M)=-BR2(J+1) 

B(N-1+J,N-1)=-2.D0*H*(BR1(N)-2.*P/R(N)*(BR8(N)-BR8(J+1) 

&))-2.D0*BR2(N) 

B(N-1+J,N+J)=BR4(J+1)+2.D0*BR4(N) 

B(N-1+J,N-1+J)=2.D0*H*BR3(J+1) 

B(N- 1 + J, N-2 + J)=-BR4( J+ 1 ) 
B(N-1+J,2*N-2)=-2.D0*H*BR3(N)-2.D0*BR4(N) 

B(N-1 +J,2*N- 1 +J)=2.DO*H*BR7(J+ 1 ) 
B(N-l+J,3*N-2)=-2.D0*H*BR7(N) 

C 

C(N- 1 +J, J+ 1 )=BR6(J+ 1 )+2.D0*BR6(N) 
C(N-1+J,J)=2.D0*H*BR5(J+1) 

C(N-1+JJ-1)=-BR6(J+1) 

C 

C AT J=N-2 
J=N-2 

A(N-1+J,J+1)=BL1(J+1)-2.D0»BL1(N) 

A(N-1+JJ-1)=-BL1(J+1) 

A(N- 1 +J,N-2)=2 . D0*BL 1 (N) 

A(N-l+J,2*N-3)=2.DO*H*BL2(J+l) 

A(N-l+J,2*N-2)=-2.D0*H*BL2(N) 

C 

B(N-1+J,N-1)=BR2(J+1)-2.D0*H*(BR1(N)-2.*P/R(N)* 

&(BR8(N)-BR8(J+1)))-2.D0*BR2(N) 

B(N-1+J,J)=2.D0*H*BR1(J+1)+2.D0*BR2(N) 

B(N-1+J,J-1)=-BR2(J+1) 

B(N-1+J,N+J)=BR4(J+1)-2.D0»H*BR3(N)-2.D0*BR4(N) 

B(N-1+J,N-1+J)=2.D0*H*BR3(J+1)+2.D0*BR4(N) 

B(N- 1 +J,N-2+J)=-BR4( J+ 1 ) 

B(N-l+J,2*N-l+J)=2.DO*H*BR7(J+l) 

B(N-l+J,3*N-2)=-2.D0*H*BR7(N) 

C 

C(N-1+J,J)=2.D0*H*BR5(J+1)+2.D0*BR6(N) 
C(N-1+J,M)=-BR6(J+1) 

(2*N-2)TH EQUATION FROM BOUNDARY CONDITION 

A(2*N-2,N-I)=H 
B(2*N-2,N-2)=SF/(R(N)*W(N)) 

B(2*N-2,N-I)=-SF/(R(N)*W(N)) 



COMPONENTS OF A,B & C FROM SIGMA II*** 
AT J=I 

A(2*N-1,I)=DL1(1) 

A(2 *N- 1 ,2 *N- 1 )=H*DL3 ( I ) 

C 

B(2*N-I,I)=DRI(I) 

B(2*N-1,N)=DR3(1) 



n o 



130 



B(2*N-1.2*N)=DR7(1) 

B(2*N-1,2*N-1)=H*DR6(1)-DR7(1) 

C 

C(2*N-1,I)=DR5(1) 



AT J=2 

A(2*N,2)=DL1(2) 
A(2*N,N)=2.D0*H*DL2(2) 
A(2*N,2*N)=2.D0*H*DL3(2) 

C 

B(2*N,2)=DR1(2) 

B(2*N,N)=2.D0*H*DR2(2) 

B(2*N,N+1)=DR3(2) 

B(2*N,2*N-1)=-DR7(2) 

B(2*N,2*N)=2.D0*H*DR6(2) 

B(2*N,2*N+1)=DR7(2) 

B(2*N,4*N-2)=2.DO*H*DR8(2) 

C 

C(2*N, l)=2.DO*H*DR4(2) 
C(2*N,2)=DR5(2) 

C 

DO 62 J=2,N-3 

A(2 *N- 1 + J, J- 1 )=-DL 1 (J+ 1 ) 

A(2*N-1+J,J+1)=DL1(J+1) 

A(2»N- 1 +J,N- 1 +J)=2.D0*H*DL2(J+ 1 ) 
A(2*N-1+J,2*N-1+J)=2.D0*H*DL3(J+1) 
C 

B(2 *N- 1 +J,M )=-DR 1 (J+ 1) 

B(2*N- 1+J, J+ 1 )=DR 1 (J+ 1 ) 

B(2 *N- 1 + J,N-2+J)=-DR3 (J+ 1 ) 
B(2*N-1+J,N-1+J)=2.D0*H*DR2(J+1) 
B(2 *N- 1 + J,N+ J)=DR3 ( J+ 1 ) 

B(2 *N- 1 +J,2 *N-2+J)=-DR7( J+ 1 ) 
B(2*N-1+J,2*N-1+J)=2.D0*H*DR6(J+1) 
B(2 *N- 1 + J,2 *N+J)=DR7(J+ 1 ) 
B(2*N-l+J,4*N-3+J)=2.D0*H*DR8(J+l) 
C 

C(2 *N- 1 + J,M )=-DR5( J+ 1) 

C(2*N-1+J,J)=2.D0*DR4(J+1) 

C(2*N-1+J,J+1)=DR5(J+1) 

62 CONTINUE 
C 

J=N-2 

A(2 *N-1 + J, J- 1 )=-DL 1 (J+ 1 ) 

A(2 *N- 1 +J, J+ 1 )=DL 1 (J+ 1 ) 

A(2 *N- 1 + J,N- 1 + J)=2 . *H*DL2(J+ 1 ) 
A(2*N-1+J,2*N-1+J)=2.D0*H*DL3(J+1) 
C 

B(2 *N- 1 +J, J- 1 )=-DR 1 (J+ 1 ) 
B(2*N-1+J,J+1)=DR1(J+1) 

B(2 *N- 1 +J.N-2+ J)=-DR3 (J+ 1 ) 
B(2*N-1+J,N-1+J)=2.D0*H*DR2(J+1) 

B(2 *N- 1 + J,N+ J)=DR3 (J+ 1 ) 

B(2 *N-1 + J,2 *N-2+ J)=-DR7(J+ 1 ) 



n n n n n 




B(2*N-1+J,2*N- 1 +J)=2.D0*H*DR6(J+ 1) 
B(2 *N- 1 +J,2*N+ J)=DR7(J+ 1 ) 
B(2*N-1+J,4*N-3+J)=2.D0*H*DR8(J+1) 
C 

C(2 *N- 1 +J, J- 1 )=-DR5(J+ 1 ) 
C(2*N-1+J,J)=2.D0*DR4(J+1) 

C 

C ATJ=N-1 

A(3*N-2,N-2)=-DLl(N) 

A(3*N-2,N-1)=DL1(N) 

A(3*N-2,2*N-2)=H*DL2(N) 

A(3*N-2,3*N-2)=H*DL3(N) 

C 

B(3*N-2,N-1)=DR1(N) 

B(3*N-2,N-2)=-DRl(N) 

B(3*N-2,2*N-2)=H*DR2(N)+DR3(N) 

B(3*N-2,2*N-3)=-DR3(N) 

B(3 *N-2,3 *N-2)=H*DR6(N)+DR7(N) 

B(3*N-2,3*N-3)=-DR7(N) 

B(3*N-2,5*N-4)=H*DR8(N) 

C 

C(3*N-2,N-2)=-DR5(N) 



*** COMPONENTS OF A,B & C FROM SIGMA33 
AT J=1 

A(3*N-1,1)=2.D0*H*EL1(2) 
A(3*N-l,N)=2.DO*H*EL2(2) 
A(3*N-l,3*N-l)=2.DO*H*EL3(2) 

C 

B(3*N-1,1)=2.D0*H*ER1(2) 

B(3*N-1,2)=ER2(2) 

B(3*N-l,N)=2.DO*H*ER3(2) 

B(3*N-1,N+1)=ER4(2) 

B(3*N-1,3*N-1)=2.D0*H*ER6(2) 

B(3*N-1,3*N)=ER7(2) 

B(3*N-l,4*N-2)=2.DO*H*ER8(2) 

C 

C(3*N-1,1)=2.D0*H*ER5(2) 

C 

DO 63 J=2,N-2 

A(3*N-2+JJ)=2.D0*H*EL1(J+1) 

A(3*N-2+J,N-1+J)=2.D0*H*EL2(J+I) 

A(3 *N-2+J,3 *N-2+ J)=2.D0*H*EL3(J+ 1 ) 

C 

B(3 *N-2+ J, J+ 1 )=ER2 (J+ 1) 

B(3*N-2+J,J)=2.D0*H*ERl(J+l) 

B(3*N-2+J,J-l)=-ER2(J+l) 

B(3*N-2+J,N+J)=ER4(J+l) 

B(3*N-2+J,N-1+J)=2.D0*H*ER3(J+1) 

B(3 *N-2+J,N-2+J)=-ER4(J+ 1 ) 

B(3 *N-2+J,3 *N-3+J)=-ER7(J+ 1 ) 

B(3 *N-2+J,3 ♦N-2+J)=2.D0’*H*ER6(J+ 1 ) 




B(3 *N-2+J,3 *N- 1 +J)=ER7(J+ 1) 
B(3*N-2+J,4*N-3+J)=2.DO*H*ER8(J+l) 

C 

C(3*N-2+J.J)=2.D0*H*ER5(J+l) 

63 CONTINUE 
C 

C ATJ=N-1 

A(4*N-3,N-1)=H*EL1(N) 

A(4*N-3,2*N-2)=H*EL2(N) 

A(4'»N-3,4*N-3)=H*EL3(N) 

C 

B(4*N-3,N-1)=H*ER1(N)+ER2(N) 

B(4*N-3,N-2)=-ER2(N) 

B(4*N-3,2*N-2)=H*ER3(N)+ER4(N) 

B(4*N-3,2*N-3)=-ER4(N) 

B(4*N-3,4*N-3)=H*ER6(N)+ER7(N) 

B(4*N-3,4*N-4)=-ER7(N) 

B(4*N-3,5*N-4)=H*ER8(N) 

C 

C 

C *♦* COMPONENTS OF A,B & c FROM TAU 22 *** 
C 

C ATJ=1 

A(4*N-2,N)=2.D0*H*FL1(2) 

A(4*N-2,4*N-2)=2.D0*H*FL2(2) 

C 

B(4*N-2,2)=FR1(2) 

B(4*N-2,N+1)=FR4(2) 

B(4*N-2,N)=2.DO*H*FR3(2) 

B(4*N-2,4*N-1)=FR6(2) 

B(4*N-2,4*N-2)=2.D0*H*FR5(2) 

C 

C(4*N-2,1)=2.D0*H*FR2(2) 

C 

DO 64 J=2,N-2 

A(4*N-3+J,N- 1 +J)=2.D0*H*FL 1 (J+ 1 ) 
A(4*N-3+J,4*N-3+J)=2.D0*H*FL2(J+l) 

C 

B(4*N-3+J,J+l)=FRl(J+l) 

B(4 *N-3+J, J- 1 )=-FR 1 (J+ 1 ) 

B(4 *N-3+J,N+ J)=FR4( J+ 1 ) 

B(4*N-3+J,N-l+J)=2.D0*H*FR3(J+l) 

B(4*N-3+J,N-2+J)=-FR4(J+l) 

B(4*N-3 +J,4 *N-2+J)=FR6( J+ 1 ) 

B(4*N-3+J,4*N-3+J)=2.D0*H*FR5(J+l) 

B(4*N-3+J,4*N-4+J)=-FR6(J+l) 

C 

C(4*N-3+J,J)=2.D0*H*FR2(J+1) 

64 CONTINUE 
C 

C ATJ=N-1 

A(5*N-4,2*N-2)=H*FL1(N) 

A(5*N-4,5*N-4)=H*FL2(N) 

C 



o o o n o 



133 



B(5*N-4,N-l)=FRI(N) 

B(5*N-4,N-2)=-FRl(N) 

B(5*N-4.2*N-2)=H*FR3(N)+FR4(N) 

B(5*N-4,2*N-3)=-FR4(N) 

B(5*N-4,5*N-4)=H*FR5(N)+FR6(N) 

B(5*N-4,5*N-5)=-FR6(N) 



COMPONENTS OF AA.BB & CC FROM RADIAL FORCE BALANCE *** 
AT J=1 

AA(1,1)=2.D0*H**2*CL1(2) 

AA(1,2)=H*CL2(2) 

AA(1,N)=2.D0*H**2*CL3(2) 

C 

BB(1,1)=2.D0*H**2*CR1(2)-4.D0*CR3(2) 

BB(1,2)=H*CR2(2)+2.D0*CR3(2) 

BB(1,N-1)=-4.*H**2*P/R(N)*CR10(2) 

BB(1,N)=2.D0*H**2*CR4(2) 

BB(1,N+1)=H*CR5(2) 

BB(1,2*N)=2.D0*H**2*CR8(2) 

BB(1,3*N-1)=2.D0*H**2*CR9(2) 

C 

CC(1,1)=2.D0*H**2*CR6(2) 

CC(1,2)=H*CR7(2) 

C 

DO 65 J=2,N-3 
AA(J,J-1)=-H*CL2(J+1) 

AA(J,J)=2.D0*H**2*CL1(J+1) 

AA(JJ+l)=H*CL2(J+l) 

AA(J,N-I+J)=2.D0*H**2*CL3(J+1) 

C 

BB(J,J-1)=2.D0*CR3(J+1)-H*CR2(J+1) 

BB(J,J)=2.D0*H**2*CR1(J+1)-4.D0*CR3(J+1) 

BB(J,J+1)=2.D0*CR3(J+1)+H*CR2(J+1) 

BB(J,N- 1 )=-4. *H* *2 *P/R(N)* CR 1 0( J+ 1 ) 

BB(J,N-2+J)=-H*CR5(J+ 1 ) 

BB(J,N-1+J)=2.D0*H**2*CR4(J+1) 

BB(J,N+J)=H*CR5(J+1) 

BB(J,2*N-1+J)=2.D0*H**2*CR8(J+1) 

BB(J,3*N-2+J)=2.D0*H**2*CR9(J+l) 

C 

CC(J,J-1)=-H*CR7(J+1) 

CC(J,J)=2.D0*H**2*CR6(J+I) 

CC(J,J+1)=H*CR7(J+1) 

65 CONTINUE 
C 

C ATJ=N-2 
J=N-2 

AA(N-2,N-3)=-H*CL2(J+I) 

AA(N-2,N-2)=2.D0*H* ♦2*CL I (J+ 1) 

AA(N-2,N-I)=H*CL2(J+I) 

AA(N-2,2*N-3)=2.D0*H**2*CL3(J+I) 

C 



134 



BB(N-2,N-3)=2.DO*CR3(J+l)-H*CR2(J+l) 

BB(N-2,N-2)=2.D0*H**2*CR1(J+1)-4.D0*CR3{J+1) 

BB(N-2,N-l)=2.DO*CR3(J+l)+H*CR2(J+l) 

&-4.*H**2*P/R(N)*CR10(J+l) 

BB(N-2,2*N-4)=-H*CR5(J+l) 

BB(N-2,2*N-3)=2.D0*H**2*CR4(J+1) 

BB(N-2,2*N-2)=H*CR5(J+I) 

BB(N-2,3*N-3)=2.DO*H**2*CR8(J+l) 

BB(N-2,4*N-4)=2.D0*H**2*CR9(J+1) 

C 

CC(N-2,N-3)=-H*CR7(J+l) 

CC(N-2,N-2)=2.D0*H**2*CR6(J+1) 

C 

C dTz EQUATION FOR (N-1) TH EQUATION **♦ 

C 

AA(N-1,2*N-2)=2.*H*(1.-F)*R(N) 
BB(N-1,N-1)=2.*H*(F*W(N)»S11(N)/SF-2.*(I.-F)*SF 
&*WP(N)AV(N)/R(N)) • 

BB(N-1,2*N-2)=2.*H*(F*R(N)*S11(N)/SF-(1.-F)*SF 

&*WP(N)AV(N)**2)-2.*(1.-F)/W(N)*SF 

BB(N-1,2*N-3)=2.*(1.-F)/W(N)*SF 

BB(N-1,3*N-2)=2.*H*F*R(N)*W(N)/SF 

CC(N-1,N-1)=-H 

CC(N-l,N-2)=-2.*(l-F)*R(N)*W(N) 

C 

C ***********♦*♦♦♦♦♦♦♦**♦**♦♦♦♦♦♦**♦♦♦♦*♦♦ INVERSE OF ICCI 
C 

CALL DLINRG(N-1,CC,N-1,CCI,N-1) 

C 

C RECURSION FORMULA WHEN THE INVERSE IS VERY ILL-C 
C CONDITIONED ( IS NOT LISTED HERE. SEE THE MAIN PROGRAM 
C FOR NEWTONIAN SINGLE-LAYER CASE.) 

C 

Q ♦♦♦*♦♦*♦*. ************^* MULTIPLYING [CCI](AAI 
C & (CCI][BB] 

DO 700 I=I,N-I 
DO 700 J=l,5*N-4 
CCIAA(U)=0. 

CCIBB(I,J)=0. 

DO710K=l,N-l 

CCIAA(I.J)=CCIAA(U)+CCI(I,K)*AA(K,J) 

710CCIBB(U)=CCIBB(I,J)+CCI(I,K)*BB(K,J) 

700 CONTINUE 

C 

Q ♦.*.***.**♦**♦♦**♦*♦♦** COMPUTE [CHCCIAA] AND [C)(CCIBB] 

C' 

DO 830 I=l,5*N-4 
DO 830 J=l,5*N-4 
CCCIAA(I,J)=0. 

CCCIBB(U)=0. 

DO 831 K=1,N-1 

CCCIAA(U)=CCCIAA(1,J)+C(I,K)*CCIAA(K,J) 

83 1 CCCIBBIL J)=CCCIBB(U)+C(I,K)*CCIBB(K, J) 



135 



830 CONTINUE 
C 

Q ***•«»*****«***,«*«,„ (a|=|A|-(CCCIAA] 
C ****.«*»**««4.«***,«,. (B1=1B1-(CCCIBB1 
C 



DO 860 I=l,5*N-4 
DO 861 J=l,5*N-4 
A(I,J)=A(I,J)-CCCIAA(I,J) 
861 B(I,J)=B(U)-CCC1BB(I,J) 
860 CONTINUE 






COMPUTE EIGENVALUE FROM DGVLRG 



CALL DGVLRG(5*N-4,B,5*N-4,A,5*N-4,ALPHA,BETA) 

DO 71 1=5 *N-4, 1,-1 

1F(BETA(I).NE.O.O) THEN 

E VAL(2 *1- 1 )=DBLE( ALPH A(2 • I- 1 ))/BET A(I) 

EVAL(2*I)=DBLE(ALPHA(2*I))/BETA(I) 

ELSE 



EVAL(2*I-1)=AMACH{2) 

EVAL(2*I)=AMACH(2) 

ENDIF 

71 CONTINUE 

WRITE(*,901) N,5*N-4 
WRITE(9,901) N,5*N-4 

901 FORMAT(//5X, 'NUMBER OF NODE= ',I3,5X,'MATRIX SIZE= ’.13 
&y5X,’FlLE; TLST) 

BP=BP*SF**2 

WRITE(*,902) TZ,BP,BUR,TR,T,RPI 
WRITE(9,902) TZ,BP,BUR,TR,T,RPI 

902 FORMAT(//lX,'DRAW FORCE= ',F9.3,5X,'B= ',F6.3/1X,'BUR= 
&F6.3,15X,'TR= ',F9.3/1X,'T= ’,F8.4,15X,’RPI= ',F10.6,/) 

WRITE(*,903) FL,P,F,D 
WRITE(9,903) FL,P,F,D 

903 FORMAT(lX,’FL= ',F5.2,13X,’PRESSURE= ',F6.1/ 

&lX,'FLOW RATIO= ',F4.2,6X,'LAMBDA= ’,F4.2) 

WR1TE(*,986) 

WRITE(9,986) 

986 FORMAT(//3X,'NO.',13X,'REAL',16X,'IMAG',/) 

K=0 

DO J=5*N-4,5*N-9,-l 
K=K+1 

EVR=EVAL(2*J-1) 

EVI=EVAL(2*J) 

WRITE(*,987) K,EVR,EVI 
WRITE(9,987) K,EVR,EVI 
END DO 

987 FORMAT(3X,I2,2X,2E20.9) 

STOP 

END 



REFERENCES 



1. Pearson, J. R. A., and Petrie, C. J. S., "The flow of a tubular film. Part I. 
Formal mathematical representation," J. Fluid Mech. . 40 . 1 (1970). 

2. Pearson, J. R. A., and Petrie, C. J. S., "The flow of a tubular film. Part II. 
Interpretation of the model and discussion of solutions," J. Fluid Mech. . 

42, 609 (1970). 

3. Pearson, J. R. A., and Petrie, C. J. S., "A fluid-mechanical analysis of the 
film-blowing process," Plast. polvm. . 38, 85 (1970). 

4. Middleman, S., " Fundamentals of Polymer Processing ." Chapter 10, 
McGraw-Hill, New York, (1977). 

5. Petrie, C. J. S., "Film Blowing, Blow Molding and Thermoforming," in 
" Computational Analysis of Polymer Processing ." Pearson, J. R. A., and 
Richardson, S. M. eds.. Applied Science Publishers, New York (1983). 

6. Petrie, C.J.S., "Memory effects in a non-uniform flow: A study of the 
behavior of a tubular film of viscoelastic fluid," Rheol. Acta . 12, 92 (1973). 

7. Petrie, C. J. S., "Mathematical modeling of heat transfer in film blowing: a 
case study," Plast. polvm.. December, 259 (1974). 

8. Petrie, C. J. S., "A Comparison of Theoretical Predictions with Published 
Experimental Measurements on the Blown Film Process," AIChE J. . 21, 

275 (1975). 

9. Han, C. D., and Park, J. Y., "Studies on Blown Film Extrusion. II. Analysis 
of the Deformation and Heat Transfer Process," J. Appl. Polvm. Sci. . 19, 
3277 (1975). 

10. Kwack, T. H., and Han, C. D., "Development of Crystalline Structure 
during Tubular Film Blowing of Low-Density Polyethylene," J. Appl. 

Polvm. Sci. . 35 . 363 (1988). 

11. Farber, R., and Dealy, J., "Strain History of the Melt in Film Blowing," 
Polvm. Eng. Sci. . 14, 435 (1974). 



136 



137 



12. Luo, X.-L., and Tanner, R.I., "A Computer Study of Film Blowing," Polvm, 
Eng. Sci. . 25,. 620 (1985). 

13. Cain, J.J., A simulation of the film blowing process, Ph.D. thesis, University 
of California, Berkeley (1987). 

14. Cain, J. J., and Denn, M. M., "Multiplicities and Instabilities in Film 
Blowing," Polvm. Eng. Sci. . 28,. 1527 (1988). 

15. Seo, Y., and Wissler, E. H., "The Effect of Extrudate Swell on Modeling 
the Film Blowing Process," Polvm. Eng. Sci. . 29,. 722 (1989). 

16. Campbell, G. A., and Cao, B., "The Interaction of Crystallinity, 
Elastoplasticity, and a Two-Phase Model on Blown Film Bubble Shape," 
Plastic Film & Sheet. . 3, 158 (1987). 

17. Cao, B., and Campbell, G. A., "Viscoplastic-Elastic Modeling of Tubular 
Blown Film Processing," AIChE J. . 36, 420 (1990). 



18. Siegmann, A., and Nir, Y., "Structure-Property Relationship in Blends of 
Linear Low- and Conventional Low-Density Polyethylene as Blown Film," 
Polvm. Eng. Sci. . 27, 1182 (1987). 

19. Montes, S. A., and Ramirez, J. L., "The Use of Melt Strength of Low 
Density Polyethylene as a Processability Measure in Film Blowing 
Extrusion," Polvm. Eng. Sci. . 26, 388 (1986). 

20. Desper, C. R., "Structure and Properties of Extruded Polyethylene Film," T 
Appl. Polvm. Sci.. 13, 169 (1969). 

21. Maddams, W. F. and Preedy, J. E., J. Appl. Polvm. Sci. . "X-ray Diffraction 
Orientations Studies on Blown Polyethylene Films, I. Preliminary 
Measurements," 22, 2721; "X-ray Diffraction Orientations Studies on Blown 
Polyethylene Films, II. Measurements on Films from a Commercial 
Blowing Unit," 22, 2739; "X-ray Diffraction Orientations Studies on Blown 
Polyethylene Films,III. High-Stress Crystalline Orientation," 22 2751 
(1978). 



Choi, K. J., Spruiell, J. E., and White, J. L., "Orientation and Morphology 
of High-Density Polyethylene Film Produced by the Tubular Blowing 
Method and its Relationship to Process Conditions," J. polvm. Sci.. polvm. 
Phvs. Ed. . 20, 27 (1982). 



22 . 



138 



23. Yeow, Y.L., "Stability of tubular film flow: a model of the film-blowing 
process," J. Fluid Mech. . 75 . 577 (1976). 

24. Denn, M. M., Petrie, C. J. S., and Avenas, P., "Mechanics of Steady 
Spinning of a Viscoelastic Liquid," AIChE J. . 21, 791 (1975). 

25. Keunings, R., Crochet, M. J., and Denn, M. M., "Profile Development in 
Continuous Drawing of Viscoelastic Liquids," Ind. Eng. Chem. Fundam. . 22, 
347 (1983). 

26. Schrenk, W. J., and AJfrey, T., "Coextruded Multilayer Polymer Films 
and Sheets," in "Polymer Blends", D. R. Paul, and S. Newman, eds.. 
Academic Press, New York (1978). 

27. Park, C.-W., "Extensional Flow of a Two-Phase Fiber," AIChE J. . 36 . 197 
(1990). 

28. Park, C.-W., "A Study on Bicomponent Two-Layer Slot Cast 
Coextrusion," Polvm. Eng. Sci. . 31,- 197 (1991). 

29. Han, C. D., and Park, J. Y., "Studies on Blown Film Extrusion. III. Bubble 
Instability," J. AddI. Polvm. Sci. . 19, 3291 (1975). 

30. Han, C. D., and Shetty, R., "Flow Instability in Tubular Film Blowing. 1. 
Experimental Study," Ind. Eng. Chem. Fundam. . 16, 49 (1977). 

31. Kanai, T., and White, J. L., "Kinematics, Dynamics and Stability of the 
Tubular Film Extrusion of Various Polyethylenes," Polvm. Eng. Sci. . 24, 

1185 (1984). 

32. Minoshima, W., and White, J. L., "Instability Phenomena In Tubular 
Film, And Melt Spinning Of Rheologycally Characterized High Density, 

Low Density And Linear Low Density Polyethylene," J. Non-Newtonian 
Fluid Mech. , 19, 275 (1986). 

33. IMSL User's Manual, Versionl.O, IMSL Inc., Houston, Texas (1987). 

34. Hyun, J. C., "Theory of Draw Resonance," AIChE. J. . 24, 418 
(1978) 

35. Denn, M. M., " On an Analysis of Draw Resonance by Hyun," 

AIChE. L . 26, 292 (1980) 

36. Tadmor, Z., and Gogos, C. G., " Principles of polymer processing " 

Chapter 6, John Wiley and Sons, NewYork (1979) 



BIOGRAPHICAL SKETCH 



Kun-Sang Yoon was born on June 25, I960, in Korea. He received B.S, in 
chemical engineering from the Seoul National University, Seoul, Korea in 1983, and 
M.S. in chemical engineering from the Korea Advanced Institute of Science and 
Technology, Seoul, Korea, in 1985. After working as a research engineer at 
Sunkyong Industries R&D Center for five years, he joined the Department of 
Chemical Engineering at the University of Florida in 1989. His career interest is in 
synthesis and processing of polymer. 



139 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 




Chang<W^n Park, Chairman 
Assistam Professor of 
Chemical Engineering 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 




Professor of 
Chemical Engineering 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 




Chemical Engineering 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 



Ranga Narayanan 
Professor of 
Chemical Eneineerine 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 

/ut. 

Wei Shyy 

Professor of Aerospace Engineering, 
Mechanics, and Engineering Science 




This dissertation was submitted to the Graduate Faculty of the College of 
Engineering and to the Graduate School and was accepted as partial fulfillment of the 
requirements for the degree of Doctor of Philosophy. 

May 1993 



A Winfred M. Phillips 

Dean, College of Engineering 




Madelyn M. Lockhart 
Dean, Graduate School