# 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