

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


MODELING OF THERMAL EFFECTS ON A SINGLE SPAN WEB HANDLING SYSTEM By JIAN LIU Bachelor of Science Oklahoma State University Stillwater, Oklahoma 1995 Submitted to the faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE May, 1999 MODELING OF THERMAL EFFECTS ON A SINGLE SPAN WEB HANDLING SYSTEM Thesis Approved: Thesis Adviser ii ACKNOWLEGEMENTS I wish to express my special appreciation to my thesis advisor Dr. D. Alan Tree for his invaluable guidance, inspiration and encouragement throughout both my undergraduate and graduate study at Oklahoma State University. Dr. Tree has been a true mentor who offers his most sincere advice in both my academic pursuit and life endeavor. I am also very grateful for the suggestions and support from the other members of the thesis committee, Dr. Martin S. High and Dr. Jan Wagner. Special thanks are due to the School of Chemical Engineering for furnishing me with the best knowledge and understanding of Chemical Engineering that the world can offer. Generous financial support in the past three years from the School of Chemical Engineering and the Oklahoma Center for Integrated Design and Manufacture is gratefully acknowledged. Finally, I would like to dedicate this thesis to my family. I am compelled to do so by their tremendous love, support and sacrifices. It is by God's grace that I was born in such a supportive and loving family. iii Chapter TABLE OF CONTENTS Page I. INTRODUCTION 1 II. BACKGROUND 4 2.1 Web Handling Models 4 2.2 Viscoelastic Effects in Web Handling 8 2.3 Temperature Dependence of Web Properties 11 2.4 Numerical Methods 19 2.5 Summary 21 III. MODEL DEVELOPEMNT 23 3.1 System Analysis 23 3.2 Assumptions 26 3.3 Governing Equations 28 3.3.1 Mass Balance 28 3.3.2 Force Balance 29 3.3.3 Energy Balance 32 3.3.4 Initial and Boundary Conditions 37 3.4 Nondimensionalization 38 3.5 An Analytical Solution for Constant Properties 40 3.6 Summary 45 IV. NUMERICAL METHODS AND TECHNIQUES 46 4.1 Selection of Numerical Techniques 46 4.1.1 Method of Lines and PDEONE 48 4.1.2 LSODE 51 4.1.3 Overall Algorithm Structure and Subroutine Arguments 53 4.2 Tension Calculation 59 iv  V. RESULTS AND DISCUSSION 62 5.1 Verification 63 5.2 Property Functions and Parameters 66 5.3 Results 70 5.3.1 Temperature Profile 71 5.3.2 Strains and Tensions 74 5.3.3 Grid Size 78 5.4 Parametric Studies 82 VI. CONCLUSIONS AND RECOMMENDATIONS 86 6.1 Conclusions 86 6.2 Recommendations 88 REFERENCES 90 APPENDIXES 97 APPENDIX AProgram Listings of the Main Program and User Defined Subroutines 98 APPENDIX BPhysical Properties of Several Web Materials 106 APPENDIX CParametric Studies of the Temperature Profile in the Thickness Direction 110 APPENDIX DEstimation of Viscous Dissipation 113 v LIST OF TABLES Table Page 1. Argwnents of the Three User Defined Subroutines 56 2. Arguments of LSODE 58 3. Test Case System Parameters and Dimensionless Groups 63 4. Fast Convergence of the Analytical Solution 64 5. Operating Parameters of a Typical System 68 D.I. Estimated Values of the Variables in the Energy Equation 114 vi UST OF FIGURES Rg~ P~ 1. Temperature Dependency of Polypropylene Properties 24 2. Thennal Strain Illustration 26 3. A Control Volume Cut from a Web Line (not to scale) 28 4. Elongation due to Defonnation and Thermal Expansion 30 5. A Web Span with Initial and Boundary Conditions 37 6. Node Point Layout of the MOL Example Problem 49 7. Layout of the Numerical Approach 54 8. Overall Algorithm Structure 55 9. Dimensionless Temperature Profile of Test Cases 65 10. Web System Dimension Summary 68 11. Isothermal Tension at Different Temperature 70 12. Dimensionless Temperature Profile of a Single Polypropylene Web Span 71 13. Temperature Profile ofa Single Polypropylene Web Span 72 14. Dimensionless Temperature Profile ofa PP Web Span over Thickness 73 15. Strains over a Range of Initial Temperatures 74 16. Tensions over a Range ofInitial Temperatures 76 17. Tensions over a Range of Ambient Temperatures 77 vii 18. The Nonisothennal Tension Prediction Error of Different Grid Sizes 80 19. The Thennal Expansion Prediction Error of Different Grid Sizes 81 20. The Effect ofPechlet Number on the Dimensionless Temperature 83 21. The Effect ofNusselt Number on the Dimensionless Temperature 84 c.l. Dimensionless Temperature Profiles at p=0.1 for Varioas Pe's III C.2. Dimensionless Temperature Profiles at p= 0.5 for Various Pe's III C.3. Dimensionless Temperature Profiles at P= I for Various Pe's 112 CA. Dimensionless Temperature Profiles at p= 1 for Various Nu's 112 vijj A b C l p D E EI F G h k L m NOMENCLATURE web crosssectional area one half ofweb thickness heat capacity at constant pressure, at constant volume constant pressure heat capacity of liquid state polymer diffusion coefficient Young's modulus elongational strain in transport direction tension shear modulus shear modulus in an uncrosslinked polymer in the rubbery plateau surface heat transfer coefficient thennal conductivity constant, mass flow rate endtoend length of a single repeat unit length of web power law coefficient ix M n N NBBrot NSGrot NPDE NEQ Nu p Pe q R Rc t T molecular weight of a repeat unit entanglement molecular weight (a) power law exponent, (b) number of grid points in transport direction number of vertices in hydrogensuppressed graph of repeat unit number of rotational degrees of freedom in the backbone of a repeat unit correction tenn for standard molar volume correlation total number of rotational degrees of freedom in the side group of a repeat unit Tg correlation parameter of twelve structure parameters number ofpartial differential equation(s) number of ordinary differential equation(s) Nusselt number pressure Pechlet number heat flux ideal gas law constant ratio of half web thickness over web length time temperature x u V, VI, V2 v w Xu ambient temperature temperature at web boundary temperature of the ith element (a) initial temperature, (b) reference state temperature glass transition temperature internal energy web velocity, first roller velocity, second roller velocity volume van der Waals volume web width thirteenth structure parameter for Tg Greek Symbols a (a) incline angle of web surfaces, (b) dimensionless thickness coefficient of volumetric thermal expansion (a) incline angle of web sides, (b) dimensionless length coefficient of linear thermal expansion zerothorder simple connectivity index zerothorder valence connectivity index firstorder simple connectivity index xi L1L E: y y or y = =(1) y 110 8 p a "C =(1) u Captions 1\ firstorder valence connectivity index solubility parameter difference operator length between two grid points in the transport direction total strain elastic strain thennal strain shear strain tensor shear strain rate tensor or convected time derivative of shear strain tensor magnitude (second invariant) ofthe shear strain rate tensor zeroshearrate viscosity time constant in the Maxwell model dimensionless temperature density nonnal stress viscous stress tensor contravariant convected time derivative of stress tensor Poisson's ratio vector indicator molar property xii ...... Subscripts o x, y, z at reference condition unless otherwise specified indicating X, y, zdirections, respectively Mathematical Symbols  'V I D gradient operator summation operator substantial derivative operator xiii CHAPTER I INTRODUCTION In the material manufacturing and processing industry, highly flexible, continuous thin sheets and strips are referred to as webs. Polymer films, paper, textiles and metal foils are all webs. Compared to large, discrete sheets of material, webs are easier to process, transport and less expensive to manufacture. Furthermore, web handling systems can be extensively automated to yield better product quality. Web handling is a widespread and essential process in the material manufacturing and processing industry. Because of the high speed of web lines, web handling poses many intriguing practical questions. However little systematic study was done before the inception of Web Handling Research Center (WHRC) at Oklahoma State University. The mission of the WHRC is to advance the knowledge base in technologies applicable to the transport and control of continuous strip materials through fundamental generic studies. Due to the increasing speed of webtransport systems, accurate tension control has been one of the major research areas at WHRC. Poor tension control may lead to wrinkling of the web, breakage or slackness, thereby limiting utilization of plant equipment, product quality, and profitability. A computerbased program, known as Web Transport System (WTS) (Shin 1991), has been developed for use in the analysis and design of openloop and closedloop control in webtransport systems. The initial 1 versions of WTS made some iC:~aUzdhai'~llmptionsfor the ~Y1Cbtmakri:lt,s~cS:lli:basjpure elastic behavior, no phase changes, no heat transfer, land no int.:rfaeiaima "IS,' tnmsl~r.J" The nonideal effec~s iIt the dmsron icontrobo.:f web lines inclJUde slj'1paf:.fb}1we~tra.M(.·"'·'Pl and a roller, web crosssectiond aIlIea change, composition chan~..;~.;rt>...:IJ.l1~',r~t'fJfechanr~ ~ and viscoelastic effects (Shin, 1991). The WHRC.has 'made a continuing effort to model the nGJ'ni dl iF efLet:9".in :v•.'.cb~~, transport systems and to incorporate these new models .u thl~ lat·~·stl,~':·rsi(','ll')llJf. WT~T;1S ;.. the first step to model nonideal effects, Guan (1995) developed a viscoelastic model (VEM) to study and predict tension due to viscoelastic behaviors. in WI6hhandling processes. In order thtisolate thte viscoelastic behavior and simplify'the pro~Lem",Guan continued to use a few of the ideal assumptioFls in Shin's work, including uniaxial tC'Iilsion:: and isotropic web materials. The results of Guan's work demonstn:.wd,the irTlplDTJtJodeJ,oL .'. viscoelasticity in w.eb:handling, as wellas,pointing out th .... ntfCossitt) fa~':oQd;yio1'hJ,:t~:r nonideal effects. This study was motivated by the need ·for it, nonisoth"rma,l moch:.:t,H't.J .:~~ address the issue of heat transfer. Heat transfer between the web and the environment is encountered in treatment and manufacturing processes such as cooling, drying, printing, and coating. For example, a frost line can be easily located in cooling sections where high temperature polymeric web materials undergo a phase change by losing heat to the environment. Since all material properties are functions of temperature, it is important to understand the thermaJo:n.:·,l effects on a web handling system in order to have accurate tension control. This thesis will concentrate on dealing with heat transfer in webtransport systems. The remainder ofthis thesis is divided in to five chapters. Chapter 11 includes tt,,,, ,!'; 2 relevant background and a literature review. Chapter III contains the development of the nonisothermal model from the constitutive equations. The systems of equations are then nondimensionalized and reduced for the purpose of general analysis. Nondimensionalized initial and boundary conditions are also presented for the governing equations. Chapter IV describes the numerical methods used to solve the resulting nonlinear partial differential equation, and develops the tension formulation. Chapter V presents the numerical simulation results and discussion. Finally, Chapter VI concludes the thesis, and comments on the direction of future studies. 3 CHAPTER II BACKGROUND Web handling is a significant sector in the modem materials processing industry, and has attracted some research interest. However, little has been done on modeling the impact of nonideal effects on tension control due to the complexity of developing rigorous models. This chapter first reviews the development ofweb handling models, which is followed by a discussion of viscoelastic effects in web handling. Section 2.3 summarizes the available methods to determine the temperature dependency of web properties. Finally, a brief discussion on applicable numerical methods for partial differential equations is given in Section 2.4. 2.1 Web Handling Models The earliest systematic study of web handling systems can be traced to the late 1950s. Campbell (1958) examined the idealized web propulsion, tension control and guidance of web materials in a free span. Campbell proposed a tension model based on the assumption that the web material obeyed Hooke's law. In an attempt to capture viscoelastic effects, the free span tension formulation was then adjusted to include an empirical internal coefficient of viscous friction. This rather simple adjustment did not 4 satisfactorily account for viscoelasticity, and no other attempt was made to capture nonideal effects. In addition, Campbell did not account for the well known "tension transfer" to the span of interest, which results from the material being in tension in preViOUS span. Assuming that web materials obeyed Hooke's law, Grenfell (1963) incorporated the tension transfer effect into his formulation and modeled the interaction of adjacent spans in multispan systems. Grenfell's model was capable ofdescribing dynamic conditions such as a step change and a sinusoidal variation in the speed of the end roller. Grenfell found that web tensions did not instantaneously respond to a change in speed. A time constant was used to measure the delay in tension response. Furthermore, Grenfell pointed out that slippage and change of moisture content in paper sheets could upset the web tension. Grenfell's work proposed feed back control strategies to solve these problems. Brandenburg (1976) reviewed tension studies of web lines, and developed a mathematical model to answer the increasing need for processing automation. Brandenburg continued to assume that the web behaved elastically and that no mass transfer occurred in and out ofthe web. The Brandenburg model successfully accounted for the adhesion and slippage of the web in the contact region with a roller due to different velocities in adjacent spans. Brandenburg's isothermal study showed that the tensile force, stress and strain were position independent in a free span. This conclusion is incorrect in processes affected by nonideal conditions such as viscoelasticity, heat transfer, and mass transfer. Furthermore, Brandenburg modeled the effect of changes in temperature both in 5  the span and in the contact region. Brandenburg introduced the thermal strain to the elastic strain to obtain the total strain. However, the temperature dependency of web properties was not included in this model. Heat transfer between the web and the environment was modeled using a thermal time constant and an average temperature at a given position. Brandenburg's model did not solve for the temperature profile within the web, which is needed in determining local properties in a rigorous mathematical model. Taguchi et al. (1985) realized that many processing problems, e.g. web breaks, wrinkles, cutoff errors and unstable runs, could be solved by adjustment of web tension. Taguchi's experimental work also showed that the tension transfer was significantly affected by the slippage between the web and the driven roller. A computer simulation, which assumed elastic behavior of the web, was developed to study the optimum design and control system for web tension. Shin (1991) developed a model which included seven primitive elements, i.e. unwinding roll, free span, driven roller, idle roller, dancer roller, nip roller and winding roll. Combinations of these seven elements can be used to describe a complex web handling operation. This idea resembles the "unit operation" concept in chemical engineering processes. Shin reviewed and tried to summarized the impact of nonideal conditions of web material, such as temperature change, moisture change and web viscoelasticity. Shin's unified model for a multispan system included the effect of slippage, temperature variation and moisture variation. Viscoelastic behavior was not considered. Shin argued that closed loop control is necessary to solve the tension transfer problem due to slippage. In a closed loop scheme, the web handling systems are treated as interconnected 6 subsystems made up of a roller and a web span. All subsystems speeds are measured and coordinated simultaneously by a master controller applying the unified model. Although complicated, the close loop control concept has been widely discussed in the field of web handling. The major drawback of the closed loop control is the difficulty in obtaining accurate online tension measurements. Shin's approach to include the effect of temperature variation ofthe web required that the web temperature be known and input to the model. This approach did not capture the essence of heat transfer within the web and between the environment and the web. Following Shin's work, Reid and Lin (1993) studied the dynamic behavior of multispan systems. They utilized simple fixedgain PID controllers in their simulation examples, and concluded that variablegain controllers were needed to produce good dynamic behavior over a broad range of startup conditions. There are other models developed for different aspects of web handling systems. Although the importance of nonideal effects has been known for decades, many models neglected them due to the complexity of introducing even one of these nonideal effects. In addition to the models already discussed in this section, there are similar elastic models in studies of multispan systems (Dunn, 1969; Soong and Li 1979; Young et a1., ]989a; Young et a1., 1989b; Parant et a1., 1992), tension control (Martin, 1973; Henderson et aI., 1979), tension measurements (Horst and Negin, ]992), tension oscillations (Veits et aI., 1981), lateral dynamics (Shelton and Reid, 1971 a; Shelton and Reid, 1971 b), winding mechanisms (Pfeiffer, 1966; Pfeiffer, 1968; Rand and Eriksson, ]973; Pfeiffer, 1977), and wrinkles (Gehbach et a1., 1989). A few viscoelastic studies can also be found in winding mechanics (Tramposch, ]967; Mukherjee, 1974; Lin and Westmann, 1989; 7  Kalker, 1991), and single span behavior of paper handling (Hauptmann and Cutshall 1977). Among these models, Tramposch (1967) and Parant (1992) incorporated temperature effects in their models, but made no attempt to find the temperature distribution within the web. Because ofthe highly temperature sensitive viscoelastic properties, Guan (1995) pointed out the importance of obtaining the temperature distribution to accurate tension modeling and control. There are also several papers dealing with the moisture and temperature effects on web materials. However, these studies were not in the area of modeling web handling systems. Morland and Lee (1960) worked on the stress analysis for linear viscoelastic materials with temperature variation. Byrd (1972) and Roisum (1992) both investigated the moisture effect on webs. LawaI (1989), Losi and Knauss (1992) and Back and Andersson (1993) studied the temperature effects on weblike materials. 2.2 Viscoelastic Effects in Web Handling When a viscoelastic material is deformed, it exhibits both elastic and viscous components. The elastic component is manifest by the recoverable strain energy. The viscous component results in an irrecoverable energy dissipation. The ratio of viscous to elastic response of a material depends on its microstructures and the processing conditions, which include temperature, pressure, geometry of deformation, type of force, prehistory and the time interval of deformation. Most web materials exhibit viscoelastic behavior. Many commercial web handling systems are designed to operate in the region where the materials exhibit 8 primarily elastic behavior to avoid complications in the tension control. However, there are many other processes or sections of processes in which the correct modeling and control in viscoelastic web materials are essential to product quality and production management. For example, the newsprint systems are heavily dependent on the viscoelastic properties of the wet paper webs under dynamic conditions. The effect of moisture is very significant to the wet paper viscoelasticity due to the weak hydrogen bonds between the crossing paper fibers (Brezinski, 1956). The environmental temperature detennines the maximum amount of moisture the air can hold, and thus, affects the relative humidity of the processing environment. The relative humidity of the environment, in turn, affects how much and how fast the moisture of the wet paper evaporates. Therefore, three nonideal effects, viscoelasticity, temperature and moisture are encountered simultaneously in the wet paper web handling systems. For other polymeric materials, the viscoelasticity will generally increase as the processing temperature increases (Ferry, 1970; Aklonis et aI., 1972; Ward, 1983). Viscoelastic models have been under development for more than a century. There are currently three types of rheological equations of state to describe the viscoelastic behavior, i.e. the Generalized Newtonian Fluids (GNF), the linear viscoelastic models and the nonlinear viscoelastic models. However, no single rheological equation of state can universally predict the viscoelastic responses of all known viscoelastic materials in various applied situations (Spriggs et aI., 1966; Tanner, 1983; Bird et aI., 1987). The GNFs are the simplest constitutive equations. The stress and strain tensors are related by a viscosity function, which is commonly chosen to be a power law relationship to account for nonNewtonian viscosity. The GNFs are useful in modeling 9 behavior under steady shear flows conditions.. The GNFs' inability to predict elastic effects limits their application in web handling systems where the elastic component is usually nontrivial. The General Linear Viscoelastic Fluids are a large group of rheological equations of state, which linearly combine Hooke's law and the Newton's law of viscosity. The idea of this type of models is to simulate the viscoelastic response using various parallel and series arrangements of a spring and a dashpot. The Maxwell model is a typical example that uses a series combination, and results in the following expression: '\ a . "( + Al "( =TJo y. = at = = (2.1 ) Th.e Voigt model, the generalized Maxwell model and the Jeffreys model are some of the other common models in the linear viscoelastic group. Guan (1995) pointed out that the linear viscoelastic models lack material objectivity when applied to systems in motion. These models predict different results for the same moving system depending on the coordinates chosen for the model analysis. The problem with material objectivity can be traced to the partial derivatives of the linear viscoelastic models. Since all web handling systems involve translation motion, Guan did not select a linear viscoelastic model. The nonlinear viscoelastic models overcome the material objectivity problem of the linear models by introducing convected derivatives to replace the partial derivatives. Guan (1995) selected the WhiteMetzner model (Bird et aI., 1987) for the study of viscoelastic effects on web handling systems: "( + TJ(Y) '[ = TJ(Y)y , = G =(1) =(1) 10 (2.2) D where (2.3) The subscript (1) denotes the upper convected derivative. The WhiteMetzner model is a nonlinear Maxwell model with a modification to include the nonNewtonian viscosity function. The WhiteMetzner model allows for stress relaxation through viscous deformation; preserves material objectivity, independence of results on coordinate system; provides the flexibility to model diverse materials such as polymer and paper; and has sufficient mathematical simplicity to allow for a solution (Guan, 1995). 2.3 Temperature Dependence ofWeb Properties The physical properties of web materials are functions of temperature. There are three ways to find the propertytemperature relations. The first, and most accurate way, is to experimentally measure the property functions of a particular material of interest. However, measuring all the needed property functions is not efficient in any preliminary design stage. The second way is to find published data on similar materials. Tremendous amounts of experiments have been done to investigate and document the physical properties of polymeric materials. Mark (1996) has organized the literature information on the physical properties of polymers. For example, the density of atactic poiypropylene as a function of temperature for various polymers is given along with the applicable temperature range and references: (2.4)  The temperature range for this correlation is 80ae  120D e. Eq. (2.4) was regressed from 11 (2.5)  data published by Rogers (1993). Although not given over a range of temperature, the thermal conductivity and heat capacity of polypropylene can also be found under other sections of Mark's work. Other properties such as the elastic modulus and the power law coefficients are not available through Mark's work. Heat capacity and thermal conductivity data can be found in Steere's work (1966) among many other published works (Harper, 1975; Hilado, 1969). The elastic modulus can be obtained from Reding (1958). This raises the problem of searching through a large amount of literature for the data in the desired form. Since the discovery of polymeric materials, many models have been developed to predict polYmer properties. As polymers with relatively simple repeat unit structure have reached their limits, more new polYmers are being synthesized for different applications. Some models are capable of making predictions even before a polYmer is synthesized. One comprehensive and very successful prediction is the group contribution techniques in the qualitative structureproperty relationships (QSPR). Van Krevelen (1976, 1990) published the classic textbook that contains a compendium of useful QSPR relationships for polymers. The group contribution techniques conceptually break down a molecule into small fragments or groups of atoms. The total value of a molar property is the summation of all the group contributions to the property of interest. The idea is illustrated by the following equation. (Total Property) :::::: I("additive" group contribution) + I("constitutive" structural terms). In Eq.(2.5), the second summation over structural terms only has to be included for certain properties where the assumption of strict addition of group contributions results in 12 p  large errors because the environment of a given group also plays an important role in the contribution of the group to that particular property. The values ofindividual group contributions can be estimated by fitting the observed values of the properties of interest from experimental data on other polymers containing the same types of fragments. Therefore, the group contribution techniques are essentially empirical and rely on statistical analysis of data. This explains the successful prediction of some properties of relative simple structure to within 1% of experimental values. In spite of its usefulness, the group contribution techniques have a few inherent limitations. The most apparent and important limitation is the need for a substantial body of experimental data for a given property in order to derive reliable values for the contributions of a large number of molecular fragments. Attempts were made to bypass the need for experimental data with theoretical models of molecular interactions. However, this approach requires extensive amounts of timeconsuming force field and quantum mechanical calculations (Hopfinger et al., 1988). Seitz (1993) developed a set of correlations to simplify the group contribution techniques. This correlation allows the prediction of a large number of "derived" physical properties of polymers from just five "fundamental" properties: molecular weight of repeat unit, endtoend length ofrepeat unit in full extension, Van der Waals volume of the repeat unit, cohesive energy, and a parameter related to the number of rotational degrees of freedom of the backbone of a polymer chain. This approach gives reasonably accurate predictions, provided the five fundamental properties can all be estimated accurately. The major drawback of Seitz's simplification is the inability to estimate the cohesive energy (Bicerano, 1993). 13 •  Bicerano (1993) suggested an alternative approach, called the topological technique. The topology is the pattern of interconnections between atoms. Since only about a dozen nonmetallic atoms (e.g. C, Hand 0) made up a large number of common macromolecules and there are only certain spatial bonding arrangements of these atoms. The development and fonnalism of "connectivity indices" can be used to mathematically predict the properties ofmacromolecules. There are significant improvements, especially in model simplicity, over other contemporary methodologies. The detailed treatment and application of the connectivity indices to different physical properties are beyond the scope of this thesis. One can find the relevant derivations in the book published by Bicerano (1993). Bicerano first obtained the four connectivity indices from experimental data. These four connectivity indices are the zerothorder simple connectivity index ox., the zerothorder valence connectivity index 0x.v , the first order simple connectivity index lX, and the first order valence connectivity index IXV. The values of these connectivity indices are tabulated for many common polymers by Bicerano. For example, the connectivity indices for polyethylene are 1.4142, 1.4142, 1.0000 and 1.000 respectively, and for polypropylene 2.2845,2.2845, 1.3938 and 1.3938. The standard properties at room temperature are then expressed as functions ofthe connectivity indices or other standard properties. The properties at any temperature are calculated from a correlation based on the standard properties. The coefficients of these correlations are obtained empirically, i.e. by curve fitting experimental data. The results of the topological approach to find the temperature correlation of the needed properties are presented below, followed by a brief discussion of a few limitations of the topology 14 technique. The temperature dependence of the molar volume of polymers are given in two sets of correlations. For polymers with Tg ~ 298 K, 1.42 Tg + 0.15 T V(T) =V(298 K) 1.42 T g + 44.7 ' (for T s Tg and Tg ~ 298 K); (2.6) 1.57 Tg +0.30 (T  Tg ) V(T) = V(298 K) , (for T > Tg~ 298 K). (2.7) 1.42 Tg + 44.7 For polymers with Tg s 298 K, (for T < Tg < 298 K); (2.8) Vet) = V(298 K) [1 + a r (298 K) (T  298)], (for T ~ Tg and Tg < 298 K), (2.9) where aT (298 K) = ( . ) , 298 + 4.23 Tg (for Tg < 298 K). (2.10) The molar volume at room temperature V(298 K) can be determined by a correlation of the four connectivity indices and a correction term NMy. V(298 K) =3.642770 Ox +9.79897 °X v 8.542819 IX + 31.693912 IXv + 0.978655 N MY' The value ofNMv is 0 for both polyethylene and polypropylene. (2.11 )  The glass transition temperature of different polymers are fitted into a correlation of four parameters, which is tabulated in Table 6.2 of Bicerano (1993). 15  N'Ig Tg =351.00+5.63 8+31.68N 23.94 x 13 . (2.12) The four parameters are the solubility parameter 8, the correlation parameter of twelve structure parameters NTg, the number of vertices in the hydrogensuppressed graph of the repeat unit N, and the thirteenth structure parameter Xl3. The tabulated values of parameters 8, NTg, N and XI3 are 15.9, 16, 2, 0.0, respectively, for polyethylene; and 17.0, 20,3, 0.0, respectively, for polypropylene. The heat capacity of solid polymers as a function oftemperature is commonly extrapolated based on the heat capacity of the polymer at a reference temperature. Bicerano used the heat capacity at room temperature, which can be predicted using a correlation of four parameters. Only the correlations for "liquid" state polymers were chosen because most polymeric materials are in the rubbery state in web handling applications. C~ (T) = C~ (298 K) [1 + 1.3 x 103 (T  298)]. C~ (298 K) =8.162061 Ox + 23.215188 Ox v + 8.477370 NBBrOl + 5.350331 N SGroI' (2.13) (2.14) where NBSro! is the number of rotational degrees of freedom in the backbone of a repeat unit, and NSGrot is the total number of rotational degrees of freedom in the side groups of a repeat unit. For polyethylene and polypropylene, the tabulated values ofNBBro! arld NSGrol are 2.0, 0 and 2.0, 1. Other than the above mentioned characteristic material properties (e.g. Tg, V, Cp), all mechanical properties of polymers (e.g. the Young's modulus) depend on the 16 > specimen properties which include crystallinity, orientation, crosslinking, thermal history and the type of force applied. This means predictions of mechanical properties are very specimen dependent. Experimental measurement is recommended to obtain reliable correlations of mechanical properties as a function of temperature. However, Bicerano's predictions are presented as a "quick and dirty" estimation for the Young's modulus E of uncrosslinked polymers under small strain in the temperature range ofthe rubbery plateau. The shear modulus ONo(T) of an uncrosslinked polymer in the rubbery plateau, is usually assumed to be determined by the physical interactions caused by entanglement between polymer chains. With the introduction of the entanglement molecular weight Me, GNO(T) can be expressed as a function of the gas constant, density, temperature and Me: G~(T) = p(T) R T. Me Me can be very roughly estimated as: N V M Me =1039.7 + 1.36411 X to23 BBrot 3 w 1m (2.15) (2.16) where NBBrot is the number of rotational degrees of freedom in the backbone of a repeat unit; M is the molecular weight of a repeat unit; Vw is the van der Waals volume; and 1111 is the endtoend length of a single repeat unit. The estimated values of these parameters are 2, 28 g/mole, 20.5 cc/mole, 2.52 x to8 em, respectively, for polyethylene and 2,44 g/mole, 30.75 cc/mole, 2.52 x to8 em, respectively, for polypropylene. For polymers in the rubbery region (T> Tg + 300 e), the Poisson's ratio, v, is very close to 0.5. The Young's modulus relationship with the shear modulus and the 17 • Poisson's ratio results in the following simple expression for uncrosslinked polymers in the rubbery region. (2.17) The thermal conductivity of amorphous polymers can be correlated from the value at Tg. (T> Tg). (2.18) lfthe thermal conductivity at any temperature is known, k(Tg) can be calculated using Eg. (2.17) and, hence, at all other temperatures. Therefore, k(298 K) is used as the reference value, and further estimation functions are needed. II [C~(298 K)] [ UR ]3 k(298 K) =5 ·10 " V(298 K) . V(298 K) , (Tg <298 K). (2.19) The molar volume and standard liquid heat capacity that appear in Eq. (2.19) are given by Eq. (2.6) to Eq. (2.9) and Eg. (2.13). The new parameter UR is called the molar Rau function and has units of cmIO /3/(sJl3mole). The estimated value of UR is 1760 for polyethylene and 2730 for polypropylene. Bicerano (1993) did not give correlations for the power law coefficients of the nonNewtonian viscosity, i.e. m and n in Eq. (2.3). Baird & CoUias (1995) have suggested the following empirical relationship between temperature and the power law coefficients. (2.20) (2.21 )  18  The reference parameters mOand nO, and the constants Band C have to be determined from experimental data. There are two limitations of the Bicerano correlations for web materials. These correlations are mainly developed for isotropic, amorphous, linear polymers. The effects of polymer crystallinity, orientation, long chain branching and crosslinking on polymer properties are not incorporated. Secondly, the connectivity indices and the correlations are obtained empirically by considering the correlation of the topological features of molecules with the properties of interest. However, the topological approach does provide some improvement over the group contribution techniques (Bicerano 1993). Polymer properties are very structure dependent. Collecting experimental data on the specific polymer of interest and curve fitting the experimental data is always the most accurate way to estimate the temperature dependence of the properties. Published data of similar specimens will provide insight and a relatively good estimate. If both of these methods are not plausible or a quick estimate is needed, caution should be exercised when using models for prediction; because no model can account for all the structural details of a particular specimen without losing generality. 2.4 Numerical Methods The study of a nonisotherrnal web transport system requires invoking the energy balance that usually consists two or more independent variables (position and time). The conduction tenus in the energy balance are second order partial derivatives. Since the properties are temperature dependent, the governing equation and the boundary 19  conditions are often highly coupled. Analytical solutions to this type of partial differential equation can usually be accomplished only with over simplification or rough approximation in which the accuracy is often lost. Considering the efficiency of modem age computing, many researchers have chosen the numerical approach. Many numerical techniques have been successfully developed to solve a wide variety of partial differential equations. Most techniques fall into the category of the finite difference methods. There are also other methods such as finite element, Monte Carlo, spectral, variational, series expansion and eigenfunction expansion methods. Each of these methods has its strength in a certain area of applications. For example, the finite element methods are popular in the field of solid mechanics and structural engineering due to its considerable freedom in setting the location of computational elements in highly irregular geometry (Strang and Fix, 1973; Burnett, 1987; Chandrupatla and Belegundu, 199 I). The spectral methods converge more rapidly than finite element methods, but only work well for very regular geometry with smooth functions (Gottlieb and Orszag, 1977; Canuto et aI., 1988; Boyd, 1989). Many specific techniques also work well with a particular type of physical problem. Patankar (1991) has developed a computer program, CONDUCT, for the analyses of heat conduction and duct flow heat transfer. In spite of these developments, the numerical solution of physically realistic, nonlinear partial differential equations is a complicated and highly problem dependent process which usually requires a researcher to undertake the difficult and time consuming task of developing a computer program. In order to eliminate the expensive and time consuming effort to solve nonlinear 20 > partial differential equations, Madsen and Sincovec (1975a, 1975b) presented a powerful and generalized software interface, called PDEONE. Based on the method of lines, PDEONE discretizes the space dimension and converts onedimensional partial different equations into a system oftime dependent ordinary differential equations, which can be reliably solved by a robust integrator. Melgaard and Sincovec (1981a, 1981 b) later developed another software interface, PDETWO, which is dedicated to solve twodimensional partial differential equations in rectangular coordinates. Both PDEONE and PDETWO require a robust integrator for initial value problems of systems of ordinary differential equations. The authors of these interfaces suggest the ODE integrator, GEARB, developed by Hindmarsh (1972,1973). In this study, an upgraded version of GEARB, LSODE (Hindmarsh, 1983) was used as the integrator for the systems of ordinary differential equations obtained from semidiscretization by PDEONE. Section 4.1 explains why these numerical methods were chosen and how the numerical routines function with each other. 2.5 Summary At the beginning ofthis study, the nonideal effects on modeling web handling systems were beginning to be addressed. Guan (1995) has completed the modeling of viscoelastic effects on multispan web handling systems. However, there were no nonisothermal models to capture the temperature distribution within the web material under nonisothermal conditions. Existing nonisothermal models assumed an average temperature variation, and dealt with the thermal effects separately from the viscoelastic 21   effects. A rigorous nonisothermal model, therefore, is needed to supplement Guan's rigorous viscoelastic model. The desired model should include the temperature dependency of web physical properties, and show how the temperature distribution affects the web tension. There are also difficulties in solving the complex energy equation, which is a highly nonlinear partial differential equation. These situations establish the background for this study. 22 a   CHAPTER III MODEL DEVELOPMENT This chapter contains the model development for temperature effects in web lines. In Section 3.1, the thermal effects on a single span web handling system are analyzed. The related assumptions for the single span system are given in Section 3.2. In Section 3.3, a set of general governing equations and the boundary conditions are developed. Section 3.4 presents the derivation of the dimensionless form of the governing equation and boundary conditions. In Section 3.5, an analytical solution with constant material properties is derived for the verification of the numerical results and parametric studies of the nonisothermal model. 3.1 System Analysis Heating and cooling in web transport systems are characterized by heat transfer between the environment and the web line by surface convection. Heat transfer between the environment and the web line causes a temperature gradient in the web material. The temperature gradient affects the tension in the web line in two ways: (1) by affecting the physical properties and (2) by inducing a thermal strain. For example, the Young's modulus is very sensitive to temperature. Since the tension is calculated using the 23  Young's modulus, temperature variation can significantly impact tension. In addition, some moderately temperature dependent properties, such as heat capacity and thermal conductivity, can also affect the process of heat transfer. Figure 1 shows the temperature dependency of some physical properties of polypropylene. The Young's modulus (E) is so sensitive to temperature that E must be plotted in log scale to fit on the same graph with the other properties. At high temperature (> 120°C), E is about one order magnitude smaller than at room temperature. 8.0·., p x 102 (kg/m 3 ) 8.0 • • • • • • • • • I I • • • • I I I • • • • • I I I • .....,. 7.0 4.0 log(E(psi)) 3.0 80 100 120 140 160 Temperature 20 40 60 2.0 +~_ __+____+_____l o Figure 1. Temperature Dependency of Polypropylene Properties* * Density, Mark (1996); Young's modulus, Reding (1958); Heat capacity and thermal conductivity, Steere (1966). 24 l The tension level is expected to increase drastically as a polypropylene web cools. The other physical properties are the density (p), the constant pressure heat capacity (Cp), and the thennal conductivity (k). The increase in Cp with temperature implies that polypropylene is more resistant to thennal changes at high temperatures. The increase in k with temperature indicates that polypropylene has greater thermal conductivity at higher temperature due to the increasing freedom of movement of polymer chains. The density is the least sensitive property among the four properties. However, the effect of a small density change can be a very significant change in the web tension level through the generation of a thermal strain. A decreasing p with temperature means polypropylene shrinks when cooled (negative thermal strain), and expands when heated (positive thennal strain). Temperature change directly induces a thermal strain in a constrained web even in the absence of stress (Shin 1991). Figure 2 illustrates the idea of thermal strain. lJndeT no constraint, the material is free to shrink or expand. In the unconstrai ned case of Figure 2, the thermal strain is Et, but there is no tension. If constrained at both ends, the material contracts, and experiences an increase in tension when the temperature decreases. The total strain in the constrained material after the temperature change is the sum of the thermal strain and the initial mechanical strain before the temperature change. If the constrained material is heated, it expands and relaxes the tension level. If the thermal expansion is greater than the initial mechanical strain, being unable to support compression, the web material sacks. 25  Unconstrained Constrained I ~~~ : I i : I f1~. Z I High Temperature Low Temperature Figure 2. Thermal Strain Illustration For a highly viscous web material, heating due to viscous dissipation can be significant during the stretching. However, since most web systems have small draw ratios (i.e. small velocity difference between adjacent rollers), the dissipation effect is considered negligible (Guan, 1995). Therefore, the tension of a nonisothermal web transport system is produced by thermal expansion/contraction in the form of thermal strain, and mechanical drawing in the form of elastic strain. 3.2 Assumptions The thermal effects of a nonisothermal system can be modeled under a set of generic operating conditions. These operating conditions are presented as the following assumptions. 26    Assumptions: 1. The thickness of the web is much smaller than the width and length. 2. The web material behaves purely elastically. 3. The web system is operating at steady state. 4. The web is under uniaxial tension. 5. The incline angles are very small. 6. No shearing deformations exist in the planes perpendicular to the machine direction. 7. No slip between the web and the rollers. 8. The web is isotropic. 9. Inertial, gravitational and air traction forces are negligible. 10. Viscous dissipation is insignificant compared to surface convection. II. The heat transfer environments above and below the web are identical. This study further assumes that the web material is purely elastic due to the complexity in dealing with both the thermal and viscoelastic effects. Forces such as inertial, gravitational, air traction, and longitudinal tension that can potentially affect the web tension have also been neglected. Finally, this study focuses on single open span systems. Under steady state conditions, multiple span systems can be treated as a combination of single span systems due to the assumption of elastic behavior. 27 a  3.3 Governing Equations The governing equations for the single span were derived from the mass balance, force balance, and the energy balance. 3.3.1 Mass Balance Figure 3 shows that the control volume cut from a web line has a wedge shape with small inclined angles as stated in Assumption 5. Q w y x I1z .I I z (Top View) i t. (\ m?_~·z (Side View) Figure 3. A Control Volume Cut from a Web Line (not to scale) The equation of continuity of any control volume system states that the rate of mass accumulation is equal to the rate of mass going in minus the rate of mass going out. Hence, the mass balance of this control volume can be written as: 28 ....... • p(z, t)A(z, t)v z(z, t)  p(z + I!.z, t)A(z + I!.z, t)vz(z + I!.z, t) I!.z =[p(z, t + ~t)A(z, t + L1t) p(z, t)A(z, t)]. 6t (3.1 ) a A (= 2bw) is the crosssectional area of the web in the xy plane. Dividing the entire Eq. (3.1) by I!.z and taking the limit as I1z and 6t approach zero gives (3.2) Under steady state conditions (Assumption 3), all variables are independent of time. The right hand side ofEq. (3.2) vanishes. A simple final form ofthe mass balance results: pAvz =km =constan t. 3.3.2 Force Balance From Assumption 4, the web system is under uniaxial tension in the machine (3.3) direction represented by the zaxis. In other words, the stress components in both the x and y directions are zero. Only the stressstrain relationship in the z direction is considered. (3.4) Consequently, the stressstraintension relationship is straightforward. Fz = Acr z = AEE e . (3.5)  For the nonisothermal case, the strain must be reexamined. Figure 4a illustrates a material that is originally in an unstrained/unstressed condition. Figure 4b shows the 29   same material under an extensional force, and Figure 4c shows the material under an extensional force after being heated. Figure 4 illustrates that the elongation results from both the elastic defonnation due to the applied force and the thermal expansion. In the nonisothermal case, the strain will be based on the total elongation regardless of the cause. Figure 4. Elongation due to Elastic Deformation and Thermal Expansion In the case of an open span in a web handling system, Y 2 YJ t ::.......:....=E1.=E c +E. (3.6) The thermal strain for an isotropic material will be pure elongation or contraction with no shear component (Crandall, 1959). Therefore, the thennal strain can be written as: yttl 0 xy = Yyz = Yxz = , 30 (3.7) (3.8) where ~l is the coefficient of linear thennal expansion, £t is the thennal strain, l is the thennal shear strain, !1T is the change of temperature. For isotropic materials, the coefficient of linear thennal expansion equals onethird the coefficient of volumetric thennal expansion (at). (3.9) By definition, the coefficient of volumetric thennal expansion is the fractional rate of change of the specific volume as a function of temperature. Since specific volume is the inverse of density, the coefficient of volumetric thennal expansion coefficient can be related to density, which is a function of temperature. a' ~ M~) ~H:) , p p where V is the specific volume, p is the density, (3.10) and a subscript p is used to denote that the partial derivative is taken at constant pressure. Thus, given the density as a function of temperature for any web material, one can easily compute the coefficient of linear thennal expansion. The density of web materials can generally be fitted to a power series from experimental data bz peT) =A + BT +CT2+... , 31 (3.11 )  with T representing the temperature. By combining Eg. (3.9) through (3.11) and carrying out the partial differentiation, the coefficient of linear thermal expansion can be expressed in terms oftemperature: c I 1 B+2CT+... ~ = 2 3A+BT+CT +... (3.12) Expanding Eg. (3.6) in the machine direction using the stressstrain relationship, Eq. (3.4), and the linear thermal expansion coefficient (~t), Eg. (3.7), gives 0 E =_z+Rt~T z E I' . Combining Eq. (3.5) and Eg. (3.13) yields F E =_z +Rt~T. Z AE I' (3.13) (3.14) Since it is difficult to measure the crosssectional area of the web, Eq. (3.3) is used to substitute for A. F pv E =_Z__Z +~t~T Z k E m The web tension can be calculated if the relevant material properties and the web (3.15) b temperature profile are known. Eq. (3.15) will be formulated to a summation form later in Chapter IV due to the discrete temperature distribution of numerical solutions. 3.3.3 Energy Balance The energy balance of the control volume depicted in Section 3.3.1 is the most essential governing equation of this modeling study. The beginning is the equation of thermal energy presented by Bird, et al. (1960). 32   ps pDD =(v.q)p(v,v)(1::Vv). Ot (3.16) IQ The left hand side ofthe equation represents the rate of gain of internal energy per unit volume. The tenns on the right hand side of the equations are from the following energy contributions in their respective order: (l) the rate of internal energy input by conduction per unit volume) (2) the reversible rate of internal energy increase per unit volume by compression, and (3) the irreversible rate of internal energy increase per unit volume by viscous dissipation. The internal energy ofEq. (3.16) needs to be modified in terms of heat capacity and temperature. Bird et a1. (1960) suggested rewriting internal energy and specific volume by applying the chain rule. ...  (aD)  (aD) [ (ap) ]   dU = ay TdV + aT V dT =  P +T aT v dV + CydT )  (av) (aY) dV = aT pdT + ap T dp . Since most web transport systems operate under atmospheric pressure, web materials are considered under constant pressure. Eq. (3.] 8) becomes (~~) = ~~. p Combining the definition of heat capacity at constant pressure and Eq. (3.19) gives   (aY) (ap)  dY(ap) C =C +T   =C +T _. . p v aT or  v dT or  p v v Multiplying dT through Eq. (3.20) and rearranging the resulting equation yields 33  (3.] 7) n.18) (3.19) (3.20) .. (3.21 ) • Substituting Eq. (3.21) into Eq. (3.17) eliminates the tenns associated with heat capacity at constant volume. (3.22) Eq. (3.22) is rewritten in the fonn of substantial derivatives. Density is then multiplied on both sides of the equation. The continuity equation of any control volume can be expressed as: Dp _ =p(V·v). Dt Thus, the substantial derivative of specific volume can be transformed using the continuity equation (Bird, 1960). DV o(X) I Dp p=p ==(V·\i). Dt Dt p Dt Eq. (3.23) can be further simplified to DO  DT p=p(V·\i)+C p. Dt p Dt (3.23) (3.24) (3.25) (3.26) Finally, substituting Eq. (3.26) into the Eq. (3.16) yields a useful form of the equation of thennal energy. pCp DT =(V .q)  (1::V'\i). Dt (3.27)  This is the equation of thermal energy for a material at constant pressure. Eq. 34 (3.27) needs to be further simplified with the assumptions listed in Section 3.2. First, the substantial derivatives and the vector operations are expanded in a rectangular coordinate system.  (Of Of Of OfJ Pc +v +v +v  P at X ax Yay z8z (3.28) The viscous dissipation terms can be neglected due to small velocity variation in any directi<;>n (Assumption 10). All the terms associated with shear stress on the right hand side of the equation disappear (Assumption 6). Moreover, the velocities in x and y direction are negligible. Under steady state condition (Assumption 3), Eq. (3.28) is reduced to (3.29) The heat transfer area available for surface convection on both the top and bottom of the web, parallel to the yz plane, are very large compared to that of the edges of the web paralleled to the xy plane in Figure 3. The dimensions ofa typical web span (see Table 5) are 5m in length, 0.3048m (l ft) in width and 5.08 x 104 m (0.02 in) in thickness. The ratio of the crosssectional area of the yz plane to the crosssectional area of the xz plane is 600: I. Because there is very little heat transfer area on the edges, the temperature gradient toward the edges is insignificant, and the heat conduction in the y direction can 35  also be neglected. Three terms are left in the energy equation representing the forced convection due to the motion of the web in the machine direction, the heat conduction toward the large surfaces and the heat conduction in the machine direction due to a temperature gradient. Rewriting the heat fluxes in terms of temperature gradients yields c (3.30) A dimensionless group called the Pechlet number represents the ratio of the heat transfer by forced convection to that by conduction. The Pechlet number can be obtained from a dimensional analysis ofEq. (3.30) if the thermal conductivity is assumed constant. The steps of nondimensionalization of the energy equation can be found in Section 3.4. (3.31 ) The heat transfer is dominated by forced convection rather than conduction if the Pechlet number is far greater than I (Baird, 1995). Since most web systems have speeds on the order of a few hundred to a few thousand feet per minute, the Pechlet number typically ranges from a few hundred to tens of thousands. Therefore, the heat conduction in the machine direction is negligible. The final governing energy equation for the web system described in Section 3.1 and 3.2 emerges as: ~ or 8 ( or) pC v = k . p z 8z Ox Ox 36 (3.32) 3.3.4 Initial and Boundary Conditions The temperature at the beginning of the web span is specified as To. T=To @ z=O. (3.33) The temperature gradient and the heat flux in the x direction are symmetric about the centerline ofthe thickness (Assumption 11). x T~ __ ._. . . ..__. z b L Figure 5. A Web Span with Initial and Boundary Conditions As shown in Figme 5, the zaxis is defined along the centerline of the thickness of the web span to obtain a simpler boundary condition. The heat flux in the x direction is zero at the center line, which gives the first boundary condition. The second boundary condition is on the top or the bottom smface of the web where the web loses heat to the environment by the means of surface convection. This boundary condition is established by combining Fourier's law of heat conduction and Newton's law of cooling. where 37 (3.34) > h is the surface heat transfer coefficient, Tb is the temperature of the web at the bOlllldary exposed to air, Ta is the ambient air temperature. The initial and boundary conditions for the halfthickness control volume are summarized as, I.C. B.C. or =0 Bx @ z=O, @x=O, @x=b. (3.35) 3.4 Nondimensionalization The governing energy equation, Eq. (3.32), and its initial and boundary conditions, Eg. (3.35), need to be converted to dimensionless forms for general analysis, computation implementation, and future applications. A set of dimensionless variables and their derivatives are defined for these purpose. x z TT a= a b' ~ = L' e= T T ; o a aa=Bx az aT a~=L' 88= (3.36) b ' To Ta A few useful partial derivatives result from these definitions. 38 • These partial derivatives are substituted into Eq.(3.22) to obtain the dimensionless governing energy equation: (3.3 7) (3.38)  where blL is the ratio of web dimension, i.e. the half web thickness over the length of the web span. If the thermal conductivity is independent of temperature, k can be moved out [ PCpv zbJ ofthe partial derivative, and the PecWet number k :would appear on the lefthand side of the equation, together with the ratio of web dimension. Section 3.5 contains the derivation of an analytical solution of Eq. (3.38) when all material properties are asswned constant. Substituting the dimensionless variables and the derivatives to Eq. (3.35) gives the initial and boundary conditions in dimensionless form. I.e. B.C. 8=1 as =0 oa as hb =8 oa k @ p=O, @a=O, @a=l. (3.39) Another dimensionless group, the Nusselt number naturally appears. 39 (3.40) The dimensionless governing equation, shown in Eq. (3.38), is a second order nonlinear partial differential equation with nonlinear boundary conditions, Eq (3.40). The numerical solution methods and techniques are presented in the next chapter. 3.5 An Analytical Solution for Constant Properties If the material properties are assumed independent of temperature, Eq. (3.38) becomes (3.41 ) where the Pechlet number (Pe) and the web dimension ratio (Rc = b/L) are constants. The initial and boundary conditions remain in the same form as in Eq. (3.39), but the Nusselt number (Nu) is now constant. The assumption of constant properties has made the governing energy equation linear. An analytical solution is, therefore, possible for the energy equation under this assumption. Since the dimensionless temperature depends on two variables, a and P, the coupled dependency can be separated ifeis represented by the product of two new variables that are functions of only one ofthe independent variables. The method is known as separation of variables. For a parabolic partial differential equation, it has been proven that the solution can be expressed in this manner (Finlayson, 1980). Let the two new variables be x =X(a), and y = y(p). (3.42) Applying the separation of variables gives 40 e(a,p) =X(a)· y(p). The partial derivatives ofe are obtained by differentiation ofEq. (3.43). (3.43) and ae dX =y aa da ' 00 dY ap =XdP , (3.44) Substituting Eq. (3.44) into Eq. (3.41) to replace the partial derivatives ofe yields dY d2X (Pe .Rc)XdP = Y da2 . The dependent variables are separated. Another new function is defined to split Eq. (3.45) into two separate equations. I dY I d2X 2 (Pe·Rc)==A, . Y dP X da2 The solution of the first resulting ordinary different equation, dY A,2 = Y dP Pe·Rc' is very straight forward. The pdependent component ofe is where CI is the integration constant. The second resulting equation from the separation of variables is (3.45) (3.46) (3.47) (3.48) (3.49) Eq. (3.49) can be easily solved by examining its characteristic equation. The solution is 41 p x = C2 sin(Aa) +C3 cOs(Aa) , (3.50) c where C2 and C3 are integration constants. The first boundary condition is applied to solve for the integration constant. @a=O, Since Y =0 is a trivial solution, Be dX aa=O=Y. da (3.51) Carrying out the remaining arithmetic of Eq. (3.52) gives Thus, the first tenn in Eq. (3.50) vanishes and S can be expressed as: ).2 S =X(a)· y(p) =C4e PeRc~ cos(Aa), (3.52) (3.53) (3.54) where C4 = C1 C3• C4 and Aare the two remaining unknown variables in Eq. (3.54). The second boundary condition is introduced to solve for A. @a=l, Be au =Nu·S. (3.55) Substituting 8 in Eq. (3.55) with Eq. (3.54) gives Further simplification of Eq. (3.56) yields Nu tan(A)  ' A . Since tan('A) is a periodic function, there is more than one solution to Eq. (3.57), or 42 (3.56) (3.57) n = 0,1,2 .... (3.58) Due to the multiple solutions of A, Eq. (3.54) is transformed into a series: The last unknown variable C4n can now be solved by applying the initial condition: (3.59) Eq. (3.59) becomes @ p=O, 8=1. (3.60) (3.61) The series form of Eq. (3.61) prevents further reduction, unless the orthogonality of the cosine function is noticed. The orthogonality of the cosine function states: for m * n, (3.62) Therefore, Eq. (3.61) is multiplied by CO~Ama), then integrated. (3.63) According to Eq. (3.62), only when m = n, the right hand side ofEq. (3.63) is not zero. Further reduction can be done in the following steps: (3.64) 43 and, (3.65) • Substituting Eq. (3.65) back to Eq. (3.59) yields the analytical solution of the governing energy equation with constant material properties. where (3.66) (3.67) This analytical solution was used to verify the numerical solution in Chapter V. Parametric studies based on the dimensionless groups can also be conducted using Eq. (3.66) and (3.77). Eq. (3.41) has been solved analytically by Carslaw and Jaeger (1959). Carslaw and Jaeger used an infinite series approach for a more general initial condition: @ ~=O, (3.68) After proving the orthogonality of X functions, Carslaw and Jaeger took advantage of the following equation in their derivation. Their result has a different appearance: form;t:. n. (3.69) (3.70) When the thickness (L) in Eq. (3.70) is replaced with the dimensionless length, 1, and 44  sin(Am) f(a) with the dimensionless initial condition, 1, the integral becomes If the Am Nusselt number is substituted using Eg. (3.67), rearrangement ofEq. (3.70) results in Eq. (3.66). This confirms the validity of the analytical solution. 3.6 Summary In this chapter, the governing mass balance, Eq. (3.3), and force balance, Eq. (3.14), were developed for the single open span presented in Section 3.1. These two governing equations were combined to give Eq. (3.15), which can be used to calculate the web tension given the temperature profile ofthe web. The governing energy equation together with the boundary conditions was then developed, simplified and presented in the dimensionless form, Eq. (3.38) and Eq. (3.39). The solution ofEq. (3.38) gives the temperature profile of the web material under nonisothermal condition. However, these equations represent a highly nonlinear partial differential equation problem that can only be solved by numerical methods. Under the assumption of constant properties, an analytical solution represented by Eq. (3.66) and Eq. (3.67) can be found for the governing energy equation. Although such solution is not practical in the calculation of web tension because the web properties are indeed affected by temperature, the analytical solution provides a basis for the parametric study of the heat transfer phenomena in the single span web system. 45 • • CHAPTER IV NUMERICAL METHODS AND TECHIQUES An analytical solution to the dimensionless governing equation, Eq. (3.38), can not be accomplished without the assumption of constant material properties. Since such an assumption results in the loss of the physical significance of this study, the challenge was to select a reliable and convenient numerical method. Section 4.1 discusses the selection of the numerical method and the implementation procedures. This is followed by the formulation of Eq. (3.15) in Section 4.2 to accommodate the discrete temperature distribution in tension calculation. 4.1 Selection of Numerical Techniques The dimensionless governing equation, Eq. (3.38), is a highly nonlinear partial differential equation due to the temperature dependency of the material properties in the coefficients. The governing equation has two independent variables that originated from the two spatial variables in the transport and thickness directions. The second order derivative represents the heat conduction in the thickness direction. The first order term describes the forced convection in the transport direction. Since the boundary condition in the conduction (thickness) direction is in the form of the normal gradients on the 46 = boundary, and the boundary condition in the convection (transport) direction is a specified value at the beginning end, the governing partial differential equation can be treated as an initial value problem with one spatial independent variable. As discussed earlier in Chapter n, there are a number of numerical techniques for the solution ofthis nonlinear initial value, partial differential equation. The formulation of these methods for any nontrivial problem is a complicated process that is highly problem dependent. One minor adjustment in the model development phase can render the previous numerical formulation useless, and require new derivation and new techniques. Web Transport System (WTS) is a complicated simulator utilizing the idea of seven primitive elements and capable of dealing with the viscoelastic effects under either transient or steady state conditions. As research on other nonideal effects progresses, WTS will need to provide numerical solutions to individual or systems of ordinary and partial differential equations. Combining the integration computation into a few modules that are reliable and versatile, will greatly benefit the development process by saving tremendous time and effort in solving these separate problems. PDEONE developed by Madsen and Sincovec (l975a, I975b) emerges as the natural choice for initial value Partial Differential Equations (PDEs) because PDEONE uses the wellstudied method of lines (MOL) to provide centered differencing in the spatial domain. When interfaced with a robust integrator for ordinary differential equations (ODEs), PDEONE becomes a powerful integrator interface for initial value PDEs with a wide variety of initial and boundary conditions. 47 • b 4.1.1 Method of Lines and PDEONE The method of lines converts an initial value PDE or systems of initial value PDEs into a sets of coupled first order ODEs by finite difference approximations ofthe spatial derivatives. The implementation procedure is best illustrated with an example. A onedimensional rod is initially at a temperature, TB. The temperature at one end of the rod is changed to TL when t > 0, and the other end is kept a TB. The temperature variation ofthe entire rod over time is of interest. This example problem can be expressed mathematically by the following governing equation, initial and boundary conditions. aT a2T =D' (4.1) at ax2 ' I.e. T(x,O) = TB , (4.2a) B.C. T(O,t)=TB , (4.2b) T(L,t) = TL . (4.2c) where T is the target temperature function; x is the spatial variable; t is time; D is a diffusion coefficient; TB and TL are the initial and boundary values of the target function. The first step is to discretize the spatial domain, i.e. define the node points. For illustration purpose, five node points are used for this example. The initial and boundary conditions of this problem and the node points are shown in Figure 6. The distance between two adjacent node points is d. The second step is to discretize the partial differential equation at each interior node point. The second partial derivative with respect to x at each node is approximated 48 c  x ~~ L ~ Ti+1 d Ti I , i T I i1  "" '... t Figure 6. Node Point Layout of the MOL Example Problem in terms of the neighboring node values. The partial derivative with respect to time becomes a total derivative at each node point. Applying finite difference approximations at any interior point, i, yields the generalized ODE for the interior points. ( i = 2,3,4). (4.3) The boundary conditions for this problem are then used for the values ofthe exterior node points. If the boundary conditions are given in the form of normal derivative with respect to x, the ODEs for the nodes adjacent to the exterior nodes have to be modified to satisfy the boundary conditions. For the example problem, the MOL discretization yields the following system of ODEs. 49 (4.4a) (4.4b)  (4.4c) (4.4d) (4.4e) & The Method of Lines has generated three coupled ODEs, Eg. (4.4b), (4.4c) and (4.4d), which can be integrated simultaneously using the Gear method software, LSODE. The ODE integration method varies depending on the stability of the problem. A brief discussion on LSODE's solution ofODE systems can be found in Section 4.1.2. It is worth noting that the initial condition of the example problem, is now transformed to the discrete format. at t =o. (4.5) For more complex initial value PDEs than the example problem, a more general implementation scheme is necessary to account for the unlimited number of mathematical structures. According to Madsen and Sincovec (l975a), PDEONE is capable of such a task for the following structure that represents a wide class of physically realistic problems. 1 8 (c aT,) 1 a( c 8TNPDE )) 7 ax x Dk,l ax '''''7 Ox X Dk,NPDE Ox ' k = 1,2, ... , NPDE, (4.6) where NPDE denotes the number of initial value partial differential equations; and c is 0, 1 or 2 depending on whether the problem is in Cartesian, cylindrical, or spherical 50 j52 coordinates, respectively. The spatial domain is in the interval of (a, b), and the time domain is (to, 00). The initial conditions are & X E [a, b], k = 1, 2, ... , NPDE. (4.7) The boundary conditions are All ofthe coefficient functions (Uk, ~k, Yk, DkJ, fk, and ~k) must be at least piecewise continuous functions of an of their respective variables. Dkj (k, j = 1, 2, ... , may be functions oft and T, but only functions ofx othelWise. PDEONE can automatically handle all three types of boundary conditions, i.e. Dirichlet (~k = 0), Neumann (Uk = 0), or mixed (Uk:l: 0, ~k :I: 0). In addition, the initial condition need not be consistent with the boundary conditions. In other words, PDEONE does not require the initial condition to satisfy the boundary condition as x approaches either boundary a or b. The derivation ofthe finite difference approximation use by PDEONE is beyond the scope of this thesis, the details can be found in the paper written by Madsen and Sincovec (l975a). 4.1.2 LSODE Most integrators can typically solve initial value ODE systems of the form i = 1, 2, ... , NEQ. (4.9) 51 >  ODE integrators are often classified as being designed to solve either stiff or nonstiff equations. If there are more than one firstorder differential equation, the possibility of a stiff set of equations arises. A problem is considered stiff when the dependent variables change in two or more very different scales of the independent variable. Most nonstiffintegrators used the Adams predictorcorrector methods of up to twentieth order (Gear, 1971). The advantage of the integrators using the Adams method is that they usually do not require the solution of large matrix problems in their solution process since a fixed point iteration is used to solve the nonlinear equations. This remains true even when these nonstiff integrators are used with the PDEONE interface (Madsen and Sincovec, 1975a). The main disadvantage of the Adams method is that the stability regions for the higher order fonnulas are small and hence can often require very small time steps to maintain stability and accuracy. Because even simple PDE problems can become stiff as the number of spatial mesh points increases, nonstiff integrators are not always adequate. Stiff integrators usually use methods modified from Newton's method for solving the nonlinear equations (Madsen and Sincovec, 1975a). The iterative algorithm of these modified Newton's methods uses the Jacobian matrix Of / aT, and therefore repeated solutions of nontrivial matrix problems are required. The most widely used stiff integrators are based on Gear's stiffly stable fonnulas (Gear, 1971), which are generalizations of the backward difference methods. Due to the limitation of the nonstiff method and the computation requirement of the stiff method, an integrator with both nonstiff and stiff capabilities is necessary to provide flexibility for different applications. LSODE (Livennore Solver for Ordinary Differential Equations) is a powerful and 52   robust ODE integrator developed by Hindmarsh (1983) as an upgraded version of the older GEAR and GEARB packages. LSODE solves stiff and nonstiff systems represented by Eq. (4.9). In the stiff case, LSODE uses the Backward Differentiation Fonnula (BDF) methods treating the Jacobian matrix as either a fun or a banded matrix, and as either usersupplied or internally approximated by difference quotients. The resulting linear system is solved by direct methods (LV factor/solve). Adams predictor corrector methods are used for nonstiff cases without using the Jacobian. matrix. LSODE is a collection of subroutines that perfonn the following tasks: initialize integrator parameters, check for errors, set time integration method coefficients, perfonn a single step of the integration and maintain a user specified time integration accuracy, select the proper time step size and order in a dynamic way throughout the problem, and set up and solve any resulting linear algebraic systems. LSODE requires the user to provide the right hand side functions (f) ofEq. (4.9) and a main program to set up integration parameters and initiate the ODE integrator. For the initial value PDE problem derived in Chapter III, the PDEONE interface subroutine is used to convert the PDE to a system of ODE. The right hand side functions can be fonned and evaluated by PDEONE. The next section explains the details of interfacing LSODE, PDEONE, the initial value PDE, and the boundary conditions. 4.1.3 Overall Algorithm Structure and Subroutine Arguments The numerical techniques discussed in the previous section work in a straight task pattern as shown by a layout of the numerical approach in Figure 7. The initial value 53 I Initial Value PDE Eq. (3.38) & (3.39) ~ PDEONE Semidiscretization (MOL) I A System of Initial Value ODEs ~ LSODE ODE Integration (Adams/BDF) Numerical Solution Figure 7. Layout of the Numerical Approach POE and the boundary conditions, described by Eq. (3.38) and Eq. (3.39), are semidiscretized in the PDEONE subroutine using the method of Jines. The result ofthe POEONE semidiscretization is a system of initial value ODEs. The package ODE integrator, LSODE, solves this system of initial value ODEs by either the Adams predictor corrector methods or the Backward Difference Formula. The overall algorithm structure works in the reverse order of the task pattern due to initiation calls and output data processing. Figure 8 illustrates the structure of the overall algorithm. The main program initializes the problem by defining physical variables, setting parameters, calling the integration driver routine. The main program 54 also sends the numerical integration results to the tension calculation routine, and processes the output. Main Program Initialize, Call DRIVER, ...... Tension Calculation OutDut Results LSODE DRIVER ODE INTEGRATOR ROUTINES LSODE by Hindmarsh D Diffusion Coefficients MOL ROUTINES PDEONE, PDETWO by Madsen et at .E Partial Differential Equations Figure 8. Overall Algorithm Structure lBNDRY Boundary Conditions  The driver routine of LSODE makes calls to different subroutines of the LSODE package. The LSODE subroutines that require information of the system of ODEs will make subsequent calls to PDEONE. For each time step in the integration process, the PDEONE is called to perfonn semidiscretization in the spatial domain. In the actual interface, PDEONE is passed as a functional argument in the calls from the main program 55  to the LSODE driver, and from the LSODE driver to the LSODE subroutines. The LSODE subroutines recognize the first argument as the right hand side functions of the ODE system, fin Eq. (4.9). PDEONE requires the initial value PDE to be defined in three subroutines, D, F and BNDRY. Table 1 lists the arguments and their array size used in the communication between PDEONE and the three user defined subroutines. The first subroutine, D, evaluates the diffusion coefficients Dkj in Eq. (4.6). For heat transfer problems, the diffusion coefficient is the thennal conductivity. Argument Array Size Related Subroutines t Scalar D, F, BNDRY NPDE Scalar D,F,BNDRY x NPDE D,F,BNDRY U NPDE D,F,BNDRY DVAL NPDEbyNPDE D UX NPDE F DUXX NPDE F FVAL NPDE F ALPHA NPDE BNDRY BETA NPDE BNDRY GAMMA NPDE BNDRY Table 1. Arguments of the Three User Defined Subroutines 56 pas The temperature dependent function of web thermal conductivity needs to be enclosed in this subroutine. The input arguments of subroutine D consist of the number of POEs (NPOE), the independent variables (t and x), and the dependent variable (U). The subroutine returns the value of the diffusion coefficients (OVAL). The second subroutine, F, evaluates the right hand side functions fk in Eq. (4.6). As can be seen from Eq. (4.6), subroutine F requires first and second order derivatives (UX and DUXX) of the dependent variable in addition to the input arguments in Subroutine D. PDEONE approximates the first and second derivatives starting from the boundary conditions. Subroutine F only has to contain the algebraic fonnula in the format specified by Eq. (4.6). .The last user defined subroutine is BNDRY for boundary conditions evaluation. Subroutine BNDRY requires the same input arguments as subroutine D, and returns the values of the coefficient functions in Eq. (4.8), Uk, Pk, and ¥k. PDEONE makes calls to these three user defined subroutines to semidiscretize the POE problem and generate a system of ODEs. Therefore, whenever LSODE and its subroutines call PDEONE, the argument of PDEONE has to contain the number of spatial points (NPTS) and the vector variable (UOOT) of the size NPDE by NPTS to evaluate the right hand side of the ODE system. In addition, the main program call to the LSODE driver routine must include 17 arguments that are needed for problem definition, error control, method selection and optional inputs. Hindmarsh (1983) gave the detailed definition and application of these LSODE arguments. Table 2 contains brief definitions of these 17 variables and their respective values for the nonisothermal problems developed in Chapter III. 57 > Argument Definition Value Used f RHS ofODE system PDEONE NEQ Number ofODE NPDE by NPTS (1 x 101) y Dependent variable Dimensionless temperature (U) t Time Time tout Point ofdesired output Every 11100 of total length itol Tolerance method flag 1 for global absolute tolerance rtol Relative tolerance 1.0 x 105 atol Absolute tolerance 1.0 x 105 itask Task flag 1 for normal computation istate Calculation status flag 1 for normal initiation iopt Optional input flag ofor no optional input rwork Real working array Uninitiated when call. lrw Declare rwork length 11132 for NEQ = 101 iwork Integer work array Uninitiated when call liw Declare iwork length 121 for NEQ = 101 Jac User supplied Jacobian Not used mf Integration method flag 22 for BDF wi internal Jacobian Table 2. Arguments ofLSODE 58  4.2 Tension Calculation The tension level within the free span can be calculated from Eq. (3.15). Since the numerical solution of the governing energy equation generates a discrete temperature profile for the web material, the remaining task is the application of Eq. (3.15) to each small segment in the transport direction (z). Applying Eq. (3.15) to the ith segment yields (4.10) By definition, elongational strain in the transport direction for the ith segment is (4.11 ) where Lo is the length of the segment before elongation; El is the length of the segment after elongation. Web velocity and tension bear no subscript, because web velocity is constant under steady state conditions for an elastic material. Tension is also constant within a span due to Assumptions 4 and 9. The constant km can be eliminated by Eq. (3.3). Combining Eq. (4.10) and Eq. (4.11) with the elimination of km results in (4.12) L If n is the total number of segments in the transport direction, then nLo is the overaU length before elongation. Summing the length of each segment after elongation gives the overall length, EI. S9  > (4.13) Dividing nLo throughout Eq. (4.13), then subtracting 1 from both sides leads to the total elongational strain of the web span, Ez. (4.14) where Ez can be calculated from the roller speed difference using the following equation derived by Shin (1991). v2  VI E z =::...._' VI (4.15) ·Finally, with the simplification ofEq. (4.15), Eq. (4.14) can be rearranged for the tension within a span. Eq. (4.16) contains a summation of the thermal strain of each discrete element. According to the elastic assumption, these elements act as tiny springs connected in a series pattern. The overall tension of the connection is evidently contributed by the deformation or strain of each tiny spring. (4.16)  According to the discretization of the numerical approach shown in Figure 6, the nonisothennal problem has two dimensions, x and z. The z dimension (transport direction) is treated as time in the integration because Eq. (3.38) has the form of an initial value PDE. The result of the numerical integration is a temperature distribution in two dimensions. However, Eq. (4.16) needs only the temperature distribution in the z dimension due to the uniaxial tension assumption. In order to remedy this problem, an 60   average temperature over the distribution in the x dimension (thickness direction) is used in Eq. (4.16) to meet the uniaxial tension assumption. This approximation does not consider the shear stress created by different degrees of thermal expansion from the centerline to the surface of the web. 61 ..  CHAPTER V RESULTS AND DISCUSSION The nonisothermal model developed in Chapter III is capable of predicting the web tension based on the analysis ofthe temperature profile within a free web span. Since polypropylene (PP) is one of the most common web materials, PP was used to demonstrate the effect of nonisothermal conditions in web handling. Section 5.1 compares the dimensionless temperature results of the analytical solution of two test cases to the numerical results of the nonisothermal model. The comparison verifies the accuracy of the temperature profile computed by the numerical methods discussed in Chapter IV. Section 5.2 contains the temperature functions of the polypropylene properties obtained from published experimental data, as well as a list of other operating parameters. The tension module of the nonisothermal model is also verified under isothermal conditions. Section 5.3 presents the numerical simulation results of a physically realistic case under nonisothermal conditions. The resulting temperature profile of an open web span is examined, followed by a study of the thermal effects on web strains and tensions. The grid size of the numerical solution is optimized for computation time with 0.1% error tolerance in the nonisothermal tension prediction. Section 5.4 concludes this chapter with parametric studies of the model based on the dimensionless groups to provide insight to the nonisothermal effects on a web span. 62   5.1 Verification The web handling industry traditionally obtains the web temperature by online measurements, and controls tension using online empirical data. As with many aspects of web handling, there is a remarkable lack of data in the open literature to compare to the results of this study. In order to verify the numerical results of the nonisothermal model, two simple test cases were devised. The test cases use the following sets of constant parameters for polypropylene. The calculated values of the dimensionless groups for these parameters are also listed in Table 3 respectively. Case #1 represents a typical operating environment of a web system. Case #2 has the same parameters except that the web velocity is reduced almost 20 times to give an extremely low Pe in web handling. System Parameters of Test Case Case #1 Case #2 Density (p, kg/mJ ) 910 910 Heat Capacity (Cp , J/kgK) 2140 2140 Velocity (vz, ftlmin) 400 20.534 Thermal Conductivity (k, J/msK) 0.172 0.172 Surface Heat Transfer Coefficient (h, W/m2K) 100 100 Length of Web Span (L, m) 5 5 Half Thickness of Web (b. m) 2.54 x 10'4 2.54 x 104 Dimensionless Groups of Test Case Case #1 Case #2 Pechlet Number (Pe) 5844 300 Nusselt Number (Nu) 0.1477 0.1477 Size Ratio (Rc) 5.08 x 105 5.08 X 10'5 Table 3. Test Case System Parameters and Dimensionless Groups 63   The most important part of the nonisothennal model is the solution ofthe governing energy equation to obtain the temperature profile, which is then used to calculate the web tension. lfthe material properties are assumed constant, Eq. (3.38) becomes Eq. (3.41). Eq. (3.66) and Eq. (3.67) express one fonn ofthe analytical solutions. where . () 2 00 sm Am ~~ e= L: e Pe·Rc CO~A a) m=O Am + sin(2Am) m , 2 4 (3.66) (3.67) This analytical solution is in the form of an infinite series that converges quickly with two to three tenns for ~ values between 0 to 0.2. For larger Pvalues, the series is accurate up to the fourth significant figure using just the fust tenn. ~ S (0) e(I) e(2) e(3) e(4) LoS 0.000 1.0234 0.0286 0.0074 0.0033 0.0019 1.0000 0.025 1.0113 0.0122 0.0003 0.0000 0.0000 0.9994 0.050 0.9994 0.0052 0.0000 0.0000 0.0000 0.9943 0.075 0.9876 0.0022 0.0000 0.0000 0.0000 0.9854 0.100 0.9760 0.0009 0.0000 0.0000 0.0000 0.9751 0.125 0.9645 0.0004 0.0000 0.0000 0.0000 0.9641 0.150 0.9531 0.0002 0.0000 0.0000 0.0000 0.9530 0.175 0.9419 0.0001 0.0000 0.0000 0.0000 0.9419 0.200 0.9308 0.0000 0.0000 0.0000 0.0000 0.9308 0.225 0.9199 0.0000 0.0000 0.0000 0.0000 0.91'99 0.250 0.9090 0.0000 0.0000 0.0000 0.0000 0.9090 0.275 0.8983 0.0000 0.0000 0.0000 0.0000 0.8983 0.300 0.8877 0.0000 0.0000 0.0000 0.0000 0.8877 Table 4. Fast Convergence of the Analytical Solution 64  Table 4 shows the fast convergence of the analytical solution for the Test Case #1. The number in parentheses after 8 denotes the mlh term of the infinite series. The last column is the sum of all terms in the series. The results ofthe analytical solution for the test cases are plotted in Figure 9, along with the numerical results of the nonisothennal model. 1.0 Case # I e 09 .0.1.) E ~ 8. 0.8 e01) f Surface r.n r.n 01) 07 . "'2 0 'r;;  Analytical Solution c 01) E 06 • NIT Model 0 Pe = 5884 Nu =0.1477 05 00 0.1 02 0.3 0.4 0.5 06 07 08 09 1.0 Dimensionless Web Length (P) 1.0 0.9 Case #2 e 0.8 .0.1..) .2 0.7 .r.o.. 8. 0.6 . E  Analytical Solution 01) f 0.5 r.n • NIT Mode.1 r.n 01) 0.4 "'2 0 Pe = 300 Nu =0.1477 'r;; 0.3 . C 01) E 0.2 is 0.1 0.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Dimensionless Web Length (P) Figure 9. Dimensionless Temperature Profile of Test Cases 65  L The results of the analytical solution are shown by the solid lines to indicate the continuous nature of the analytical solution. The results of the numerical nonisothermaI model are represented by the square dots. As can be seen from Figure 9, the dimensionless temperature results obtained by both methods coincide with each other. For Case # I, the average relative difference between the numerical results and the analytical results is less then 0.01%. For the extreme Case #2, the average difference of the results is about 0.62% due to the almost zero dimensionless temperature for ~ greater than 0.6. In both cases, the dimensionless temperature results, on a scale from ato I, agree to the fourth significant figure, which is the specified error tolerance level of the numerical integration. This agreement indicates that the nonisothermal model can solve the governing energy equation accurately, and has been properly programmed for the test cases with constant properties. For a realistic heat transfer problem, the nonisothermal model incorporates the property functions in its numerical solution. The next section contains these property functions. 5.2 Property Functions and Parameters Four material properties that are functions oftemperature need to be supplied to the nonisothermal model before any analysis is possible. They are density, thermal conductivity, constant pressure heat capacity, and the Young's modulus. Predictive equations to estimate the temperature dependency of these properties can be found in Section 2.3. However, first hand experimental data and literature data are preferable. For 66  (5.1 )  L the purpose oftbis study, temperature functions of these four polypropylene properties are obtained by curve fitting the literature data. The unit of temperature is °C in all the following functions. The density data can be found in Mark (1996), and the function is P(kgl m3 ) = 8481.90x 1O2T3.05x 1O3T2 . The data for constant pressure heat capacity per unit volume was retrieved from Steere (1966), and the function is Cp(Call cm3 . K) =0.3397 + 1.028 x 1O3T + 1.164 x 1O5T2 . (5.2) Also retrieved from Steere (1966), the thennal conductivity function is k((all s· cm . K) =4.197 x 104 + 1.082 X 106T  1.045 x 108T2 + 1.090 x 1010T3 . (5.3) The Young's modulus function was regressed from stiffness data published by Reding (1958), and is E(psi) =93460  2385T + 34.86T2  0.2506T3 + 6.496 x 104 T4 . (5.4) The temperature function for the linear thennal expansion coefficient, pl, can be obtained from the density function using Eq. (3.12). Other operating parameters needed for the nonisothennal analysis include the dimensions of the web span, the roller velocities, overall surface heat transfer coefficient, initial and ambient temperatures. The operating parameters of a typical web span system are listed in Table 5. The surface heat transfer coefficient (h) was assumed a typical constant value of 100 W1m2K (Baird, 1995). The surface heat transfer coefficient is not required to be a constant by the numerical methods. As long as h is a continuous function of temperature 67 .......  Operating Parameters Value Unit Length of Web Span (L) 5 m Half Thickness of Web (b) 2.54 x 104 m Width of Web (w) 12 In Initial Web Crosssectional Area (Ao) 0.24 in2 First Roller Velocity (VI) 400 ftfmin Second Roller Velocity (V2) 401 ftfmin Initial Temperature (To) 120 °C Ambient Temperature (Ta) 25 °C Surface Heat Transfer Coefficient (h) 100 W/m2K Table 5. Operating Parameters of a Typical System VI = 400 [pm V2 = 401 fpm z ~ 11..4.t.  L =5 m.~I Side View x W=12in z .................._.__ ._.__ .  ..~ z y Top View Enlarged Side View L Figure 10. Web System Dimension Summary 68 or location, the user only needs to supply an additional subroutine to evaluate h in the same manner that the thermal conductivity (k) is evaluated for the boundary condition. Figure 10 summarizes the dimensions of the example web system. The system is a single span web with two driven rollers whose velocities are 400 ftlmin and 401 ftlmin, respectively. The initial web temperature (To) immediately after the contact zone ofthe first roller is at 120°C in all simulation unless noted otherwise. Similarly, the ambient temperature (Ta) is at 25°C unless noted otherwise. The polypropylene property functions in the temperature range from 25°C to 120°C are described by Eq. (5.1) through Eq. (5.4). The effect of the nonisothermal condition on web tension was the ultimate goal of this study. Since strain is directly related to tension, both strains and tensions under isothermal and nonisothermal conditions need to be examined. Under isothermal conditions, the web system does not undergo a temperature change, there is no thermal strain, and all physical properties remain constant. The only strain is the mechanical strain due to the different roller speeds. According to Eq. (4.15), the mechanical strain is 2.5 x 10'3 for the system defined by Figure 10 and Table 5. Consequently, Eq. (4.16) can be further simplified to (5.5) L When the ambient temperature (Ta) and the initial web temperature (To) are set to equal to each other, the nonisothermal model can be used to generate isothermal tension results that can be compared to Eq. (5.5). The results of the nonisothermal model integration and the results from Eq. (5.5) were plotted in Figure 11. The agreement of these results verifies that the tension module of the nonisotherrnal model has been 69 accurately programmed. Both Eg. (5.5) and Figure 11 show that the isothennal tension is only affected by the temperature via the Young's modulus. Figure 11 also demonstrates the sensitivity of web tension on temperature. 35 30 ~ 25 gc0 'Vi 20 cv r '"E 15 . v .c (5 ~ 10 5 0 20 30 • NIT Model  Equation (5.5) 40 50 60 70 80 90 100 Web Temperature/Ambient Temperature (0C) 110 120 Figure 11. Isothennal Tension at Different Temperature 5.3 Results Numerical simulations have been conducted for the system depicted by the property functions of polypropylene and the operating parameters presented in Section 5.2. Section 5.3.1 presents the temperature profile results. Section 5.3.2 contains the results and analyses on the tensions and strains at different initial and ambient temperatures. Section 5.3.3 examines the effect of grid size on the precision of results. 70 . 5.3.1 Temperature Profile Figure 12 is a plot of the dimensionless temperature profile for the nonisothermal case. The trend of the dimensionless temperature profile is similar to that of Test Case #1 shown in Figure 9. However, the dimensionless temperatures in Figure 12 decrease more gradually due to the incorporation of variable physical properties. At the surface end point of the web, the dimensionless temperature was 0.5928 in Test Case #1 with constant properties. The dimensionless temperature at the same location was 0.6508 in the nonisothermal model with variable properties. ''Q"") co "V; I:: Q) E i5 0.7 o a =0.0 " a = 0.2 )( a = 0.4 x a =0.6 • a = 0.8 + a = 1.0 0.\ 02 OJ 04 05 06 07 0.8 09 10 06 '~~~~~~~l 00 Dimensionless Web Length (P) Figure 12. Dimensionless Temperature Profile of a Single Polypropylene Web Span The difference may be explained by arguing that the material constant properties in the isothermal case tended to favor conduction over convection, hence, the web temperature profile was steeper in Test Case #1. Figure I in Chapter III shows that the 71  density of polypropylene increases slightly as temperature decreases while the heat capacity and thermal conductivity both decrease. Eq. (3.31) shows that the Pechlet number contains all three properties. (3.31 ) When the web cools, the decrease of k and the slow increase of p compete with the decrease of Cpo For polypropylene, the combined effect is a higher Pechlet nwnber indicating stronger convection at lower temperature in the nonisothermal case. Therefore, less heat is lost to the environment and the temperature at the end is higher in Figure 12. oX=o A X = b/5 )( X = 2b/5 ]I X = 3b/5 • X = 4b/5 110 ,, u t.... d..o.) a 100 .'".. do) 0. E do) f< 90 +X=b 80 '.~ o Web Length (m) Figure 13. Temperature Profile ofa Single Polypropylene Web Span The dimensionless temperature was converted back to actual temperature in 72 1.  r Figure 13. As can be seen in both Figures 12 and 13, the temperature profile over the length of web span appears almost completely linear, indicating a high convection to conduction ratio. This can be attributed to the high transport speed of the web leading to a relatively high Pechlet number. A high b/L ratio (Rc) will have similar effects, but the dimensions of the web span are usually not controllable operating parameters. Q • • • • • • • • • • • • • • • • • • • .. .. .. .. .. .. .. • • • .. .. .. .. .. .. .. • • x x x x x x x x x '" '" '" x x x '" '" J< J< + + + + + + + + + + + + + + + + + + + 0 c c c c c c c c c c c c c 0 c c c c 0 0 0 0 0 0 0 0 0 a 0 a 0 0 0 0 0 0 a • ~ = 0,0 .. ~ =0.2 J< ~ = 0.4 + ~ = 0.6 c ~ = 0.8 a ~ =1.0 l.l ~ 1.0 ~ .l.U.. ::I <..;.j. 0.9 lU C. E ~ 0.8 til til '!U C0 'c;; 0.7 l: !U E 6 0.6 0.5 0.0 0.2 0.4 0.6 0.8 1.0  Dimensionless Web Thickness (a) Figure 14. Dimensionless Temperature Profile of a PP Web Span over Thickness Figure 14 is a plot ofthe dimensionless temperature profile in the thickn.ess direction for the nonisothennal case. The temperature gradient over the thickness tends to increase gradually from the center to the surface of the web. This pattern is typical when Newton's cooling occurs at the boundary. Figure 14 also shows that the temperature gradient over the thickness of the web is fairly small because the actual measurement of web thickness is in the range of 0.01 in. Parametric studies of the 73    dimensionless temperature profile in the thickness direction, with the Pechlet and Nusselt numbers, can be fOWld in Appendix C. 5.3.2 Strains and Tensions Eq. (3.6) stays that the total strain (Ez) is the sum of the elastic strain (E) and the thennal strain (E1 ). The elastic strain is only caused by the tension in the web material. Since web material does not support compression, the elastic strain is always positive. On the other hand, because the linear thennal expansion coefficient is always positive, the thennal strain is negative if the web is cooling. Figure 15 is a plot of the strains over a range of initial temperature for the system defined by Figure lO and Table 5. 0008 T. = 25°C 0.006 • • Elastic Strain (E e ) • • 0.004 • • • •• • • •• • • • • • • • c:: 0.002 .(;j Total Strain (Ez) due to 1:: Mechancial Drawing r:rJ 0000 AA A • • • • • 0002 Thermal Strain (EI ) • • • ·0004 . 0.006 20 ]0 40 50 60 70 80 90 100 110 120 rnitial Web Tern peratu re, To (0C) Figure 15. Strains over a Range ofInitial Temperatures 74  As can be seen in Figure 15, the thermal strain becomes more negative if the initial temperature of the web is higher. As the difference between the initial temperature and the ambient temperature becomes greater, the elastic strain (ge), which is the quantity (Ez  f;t), becomes greater due to the decreasing (increasingly negative) value of gt. The web tension calculated from the elastic strain becomes larger also. The greater web tension is the result of the thermal contraction. From a mathematical perspective, since the total strain is fixed by the constant roller speeds, Eq. (3.6) requires a larger elastic strain to compensate for the decrease of thermal strain. In Figure 15, the total strain due to mechanical drawing is a constant value of 2.5 x 10'3 , which is the same as in the isothermal Test Case #1. Figure 15 also demonstrates the opposite directions of the elastic strain and thermal strain as the initial web temperature increases. The tension in the open span is also the sum of two components: tension due to mechanical deformation and the tension due to thermal shrinkage. For the web system defined by Figure 10 and Table 5, the nonisothermal model has performed the numerical integration over the initial temperatures (To) from 25°C to 120°C. The resulting tensions are plotted in Figure 16. When To = 25°C, the tension due to thermal strain is zero because this initial temperature is the same as the ambient temperature and there is no heat transfer. The tension due to mechanical strain and the overall nonisothermal tension have the same value of 31.15 lbf. As To increases, the tension due to thermal strain first increases, and reaches its maximum at about To = 110°C, then starts to decrease. The first increasing behavior results from the increasing thermal strain due to higher initial web temperature. 75  ...... T. = 25°C I •• • • • • • • .. • • • .. .. .. .. .. .. • • • • • • • • • Nonisothermal Tension • • Tension due to Thermal Strain • • • • • Elastic Tension  35 30 25 ~ :9 20 '" t:: 0 'Vi vt:: 15 flO o 20 30 40 50 60 70 80 90 100 110 120  Figure 16. Tensions over a Range ofInitial Temperatures Figure 16 also shows that the thennal strain becomes more negative. The negative sign of the thennal strain indicates the web is shrinking, At about To = 11 DOC, the Young's modulus has decreased to a point to offset the effect of thennal shrinkage. The tension due to thermal strain, which is a product of the Young's modulus and the thennal strain, peaks out and turns downward. The rapid decreasing of the Young's modulus with temperature is also reflected by the fast decreasing trend ofthe tension due to mechanical strain. The nonisothennal tension, which is the sum of the tensions due to mechanical strain and thennal strain, first decreases with To, then starts to increase slowly at about 50°C, finally turns to decrease at about 9DoC. This pattern indicates the competing effects of the thennal strain and the Young's modulus under the nonisothennal condition. In this case study, the thermal strain increases and the Young's modulus decreases with 76  higher initial temperature. The two competing effects tend to lessen each other's impact on the nonisothermal tension at the mid temperature range ofTo, i.e. 40°C  100°C. The nonisothermal tension is expected to decrease drastically if the initial temperature is beyond 110°C, because the Young's modulus is so small that both the tensions due to mechanical strain and due to thermal strain are declining. 20 25 30 ,15 40 To = 120°C • • • • • · • • • • • • • Nonisothermal Tension .  • Tension due to Thermal Strain • Elastic Tension o 15 10 30 § 15 .Vj ~ 20 25 Figure 17. Tensions over a Range of Ambient Temperatures The effect of the ambient temperature, Ta, on the tensions has also been studied. Under the condition defined by Figure 10 and Table 5, the simulation tension results were plotted against ambient temperatures in Figure 17. The nonisothermal tension and its component tensions are less sensitive to Ta than to the initial web temperature (To), because the change of Ta is rather small when compared to the already large difference between To (= 120°C) and Ta. In other words, the change in the driving force of heat  77   transfer is small. The decreasing trends ofthe nonisothermal tension and both the components with increasing Ta can be interpreted in terms ofthe two most significant effects of the nonisothermal condition. When the level of the web temperature profile rises slightly with higher Ta., the Young's modulus of the web (E) becomes slightly smaner. Therefore, Figure17 shows that the tension due to mechanical strain decreases slowly with Ta. On the other hand, because the initial temperature difference, (To  Ta), becomes smaller with higher Ta, the temperature profile becomes flatter. Since the thermal strain indicates how far the web temperature decreases from the initial temperature, the flatter the temperature profile, the smaller the thermal strain. Consequently, as a product ofE and Et,the tension due to thermal strain decreases with increasing Ta. Instead of competing, E and Et reinforce each other when the ambient temperature rises. In Figure 17, the nonisothermaJ tension shows a faster decrease because it is the sum of the decreasing tensions due to both mechanical strain and thermal strain. 5.3.3 Grid Size In this section, the grid size used in the numerical integration and the tension calculation is examined to achieve an acceptable error tolerance. Since different web systems will have different physical properties and different error requirement, the web system and the error tolerance must be defined before adjusting the grid size. For the purpose of this study, the web system for the grid size analysis is again described by Figure 10 and Table 5. 78 The grid size (or mesh size) in this discussion is in the fonn ofNPTS x NSTEP. NPTS is the number of points over one half ofthe web thickness because the numerical integration is only perfonned from the center to the surface of the web in the thickness direction. NSTEP is the number of grip points in the transport direction. Since the transport direction is considered as time by the integration routine, LSODE, NSTEP is originated from the number of time steps. However, the actual integration steps are internally generated by PDEONE to match the specified error tolerance, which was set at an absolute value of 1.0 x 10.5 to give 4 significant figures accuracy for the dimensionless temperature. Therefore, NSTEP is only used to control the points where the integration routine will output the dimensionless temperature. The dimensionless temperature at the surface end point converged to 0.6508 for 20 meshes, i.e. NPTS = 21. NSTEP was only set at 2 for a single output point. The temperature accuracy was controlled internally by the integration algorithm. This demonstrates that the PDEONELSODE numerical integration approach is a stable and converges rapidly. However, since the tension calculation involves the summation ofthe properties at each grid point, the mesh size needs to be large enough to account for the sensitive temperature function of the Young's modulus. The nonisothermal tension of the web system defined by Figure 10 and Table 5 was solved repeatedly for different mesh sizes. For To = l20De and Ta = 25°C, the nonisothennal tension converged to 24.64 lbf. The mesh size for this value was found at 501 x 5001. Further increase of the mesh size shows no change in the nonisothennal tension up to the fourth significant figure. Therefore, 24.64 lbf is used as the reference tension to obtain the percent error of the 79  other mesh sizes. The percent error is rather a precision error that an accuracy error. F 24.64 %Error = x100% 24.64 (5.6) 06 c .2 U ~ 0.4 ... 0.. c0 'r;; 0.2 C Cl.l f< lil E 0.0 Cl.l 5 0en ·c 0.2 0Z .5... 0 t: 0.4 W '#. • • • • • • 0 • 0 " " " " " " 0 • • • •. 0 0 Q • • NPTS = II • o NPTS = 21 .. .. NPTS = 51 o NPTS = 101 + NPTS = 201 o NPTS = SOl 0.6 100 1000 10000 100000  Number of Transport Direction Grid Points (NSTEP) Figure 18. The Nonisotherrnal Tension Prediction Error of Different Grid Sizes Figure 18 is a plot of the percent error in the nonisotherrnal tension for different mesh sizes. From Figure 18, less than 0.1% error in the nonisotherrnal tension prediction can be achieved with a grid size of 51 x 501 or larger. If 0.1% error in precision is required, 51 x 501 is the recommended mesh size. The series ofNPTS = 101 (represented by the circles), and the series ofNPTS = 201 (represented by the crosses) coincide almost everywhere, except at NSTEP = SOland 2001. This slight fluctuation of the errors indicates how sensitive the tension is to the Young's modulus and the thermal strain. 80 ......  Since the nonisothermal tension is directly affected by the thermal strain, the percent error in thermal strain prediction can be used as another precision index. The reference strain is 4.053 x 103 , which is again obtained with a mesh size of 50 I x 5001. The definition of the percent error in thermal strain is similar to Eq. (5.6). • NPTS = II c NPTS = 21 '" NPTS = 51 o NPTS = 101 . + NPTS = 201 ., NPTS = 501 c • ~ • D • • • • '" D D ~ '" D D 0 0 ~ •'" 0'" 0'" 0'" 14 C I 2 .9 '0 :a~... 1.0 0.. c '(..;.i (j) 0.8 «i E ~ 0.6 .c: f ..!.::. 04 0 t: w ~ 0.2 0 0.0 סס 100 1000 10000 10 oo  Number of Transport Direction Grid Points (NSTEP) Figure 19. The Thermal Expansion Prediction Error of Different Grid Sizes Figure 19 shows that more grid points are necessary to keep the error less than 0.1 %. Consequently, the number of the grid points in the form ofNPTS x NSTEP is recommended to be at least 101 x 1001. The increased mesh size requirement for the same precision in thermal strain can be attributed to the competing effects of thermal strain and Young's modulus. It was concluded earlier from Figure 14 that these competing effects decrease each other's impact on the nonisothermal tension. 81 ......  5.4 Parametric Studies The nonisothermal model contains more than a dozen operating parameters and property functions. Individual parametric study only offers disconnected information on the general behavior ofthe nonisothermal model. Dimensionless groups of these operating parameters can better illustrate the effects of the parameters under nonisothermal conditions. During the derivation of the governing energy equation, two dimensionless groups, Pe and Nu, were encountered.  and hb Nu=k . (3.31) (3.40) b The Pechlet number (Pe) represents the ratio of the heat transfer by forced convection to conduction. In the nonisothermal model, the forced convection refers to the energy transported into and out of the span via the movement of the web. A portion of the heat brought in the web span control volume by forced convection is transferred to the web surface by conduction and lost to the environment by surface convection. The remaining portion of the heat is transported out of the control volume by the forced convection. According to Eq. (3.31), Pe is the gross representation of density, heat capacity, thermal conductivity, web velocity and thickness. The second dimensionless group is the Nusselt number (Nu) which shows the dimensionless temperature gradient averaged over the heat transfer surface. This definition is the same as the common definition ofNu, which is the ratio of surface convection to conduction. Both definitions can be visualized by combining Eq. (3.10) 82 ....  and the boundary condition on the surface as depicted by Eq. (3.39) (5.7) Similarly, the Nusselt number is the gross representation of surface heat transfer coefficient, web thickness and thennal conductivity. Since Pe and the Nu are dimensionless quantities, the dimensionless temperature is a natural choice for the parametric study. Figure 20 shows how the dimensionless temperature at the web center (ex = 0) corresponds to different Pechlet numbers. The Nusselt number used to obtain the data points on this plot is 0.1477, which was also used in Test Case #1. 09 10 Pe=8000 0.7 08 a = 0, Nu =0.1477 04 0.5 0.6 + + + 0.3 Pe=IOOO 0.1 0.2 Pe=300 Pe=50000 +•J•o.•. • I I ·li·S· · · · · · · · · · · · · · · · · · · · · · · · · •••••• .. •• Itt·.. ~ ... . + .. .. .. • • • : : • • • • • • l)e= 10000 + .. • • & .. • • • • I & ... • • • • A .. ... ..... ••••••• + • .. • • • • .. .. .. ... ... & • • + .. .. . Pe=5844· • • • .. ...... ... ......~ + ... & ... ;' • • • • + Pe=3000" .......... + ....... + ..... ... .to .... + + ... ... .. ... 1.0 CD '' eu 08 l aos l eu Co E 06 eu fo Cfl ~ t: 04 'r;; Ceu E 6 02 00 0.0 Dimensionless Web Length (13) Figure 20. The Effect ofPechlet Number on the Dimensionless Temperature A very large Pechlet number would correspond the web moving through the open  83 ......  span so fast that there is little or no time for conduction to the surface. The web temperature change is very small over the length for large Pechlet numbers due to the slow conduction. Figure 20 confIrms that the dimensionless temperature at the web center decreases slowly, and the profile is almost completely linear for Pe > 3000. For small Pechlet numbers, conduction is rapid compared to transport through the span. As shown in Figure 20, the smaller the Pechlet number, the faster the initial drop in dimensionless temperature. ln addition, the rate at which the dimensionless temperature decreases, slows at large dimensionless web length (P) because the conduction heat transfer driving force (T(b=O)  T(b=l) ) is reduced. Figure 21 shows how the dimensionless temperature profile at the web center (a = 0) is affected by the Nusselt number. The Pechlet number used to obtain the data points on this plot is 5844, which was also used in Test Case #1. 08 09 10 a = 0, Pe =5844 os 0.6 07 c c cacooooQCCnn~nnnn • • • • •• • • • • • • • • • • • • • • • • Nu=1 c c o • 0 0 0 • 0 • 0 0 • 0 c 0 • 00 0 • 0 • 0 0 • c • c Nu=O.OI 1 0 ~ I I I I I • I : : : • • • • • • • • • • • • • • • • • • • • • • • • • • • • o : 0 0 t t ~ •x •• •: : •••• •A • • • • • • • • • • • • • • • Nu=O.05 • 0 • • • z • • A A • A • • • • • • • • • 0. 0 0 •••• "'''' •••••••• Nu=O.1 • 0 0 • Z K * • * • • A • A • • • •• 0 00 Nu=O.1477 ••••••• • J( :I: K 02 0 0 Nu=5 c c c c c 0.0 0.0 01 0.2 OJ 04 e~ 0.8 ~ ll) 0 E 06 ~ en en ~ c: 04 'v; c: ll) E is Dimensionless Web Length (~) Figure 21. The Effect of Nusselt Number on the Dimensionless Temperature 84  ......  l Since the Nusselt nwnber is the ratio of surface convection to conduction, a very small Nusselt number represents that the convection on the surface is so slow that the web is almost insulated. The web temperature change is very small due to the slow convection. The result of a very small Nusselt number is similar to a very large Pechlet number in that little heat transfer occurs between the environment and the web. Figure 21 shows that the Nusselt number affects the dimensionless temperature in much the same pattern as the Pechlet number, except in an opposite manner. For a large Nusselt number, i.e. Nu> 0.5, the strong surface convection over conduction results in nonlinear temperature profiles. The quick initial drop in temperature resembles the effect of a Pechlet number less than 3000. The temperature profile analysis with the Pechlet number and the Nusselt number reveals that the heat transfer process in the web system involves three steps, forced convection, conduction and surface convection. Other than the temperature differences, each step is affected by some operating parameters and physical properties. Since the physical properties are in turn affected by the temperature, it is necessary to have a nonisothermal model capable to deal with this coupling relationship. 85 ......   CHAPTER VI CONCLUSIONS AND RECOMMENDATIONS 6.1 Conclusions This study was successful in modeling the effects of nonisothermal conditions on a single web span. The development of the nonisothennal model and numerical simulation provided better understanding of the heat transfer between the web and the environment. The significance of this work rests on the ability to predict the web temperature profile, hence, the local physical properties critical to tension control. Web handling systems that are under nonisothennal condition can use the nonisothermal model to accurately analyze the thermal effects, which the currently existing isothermal models and previous attempts can not capture. The major conclusions are summarized below. 1. A nonisothennal model was developed to account for the thermal effects in web handling systems whose temperature is different from the temperature of the environment. The temperature dependence of the needed material properties were incorporated into the nonisothennal model, which was verified in a test case using constant material properties. The resulting partial differential equation was solved using widely accepted and robust numerical routines, PDEONE and LSODE. 86 r  2. A modeling framework was established to accurately evaluate the temperature distribution within the web. Since the energy equation was solved independently, the model can be adjusted to solve problems with other nonideal effects in conjunction with other models, such as Guan's Viscoelastic Model (VEM). 3. Web tension was found to be significantly affected under nonisothennal conditions in two ways, by inducing a thermal strain and altering the physical properties. The larger the temperature difference between the web and the environment, the greater the predicted isothermal tension deviation from the nonisothermal tension. 4. Among the physical properties encountered in this study, the Young's modulus was the most sensitive to temperature. Since tension is directly proportional to the Young's modulus in an elastic model, accurate determination of the Young's modulus throughout the web is imperative. 5. The mesh size in the numerical simulation is as important as an accurate temperature function of the Young's modulus. Significant error in tension and strain prediction can arise if the mesh size is not large enough. Trial runs with different mesh sizes are recommended to meet specific precision requirement. 6. In a cooling process, the thermal strain is negative and the web is shrinking. The larger overall strain makes the nonisothermal tension larger than the mechanical tension calculated from the velocity difference of the rollers. In a heating process, the thermal strain is positive and the web is expanding. In order to avoid sagging, the mechanical strain has to be greater than the thermal strain to retain positive tension. 7. Two dimensionless groups, Pe and Nu, can be used to study the effect of operating parameters in the nonisothennal model over the temperature distribution. The 87 .... thermal effects on the web not only depend on the material properties and the temperature difference, but can also be controlled by the operating parameters. 6.2 Recommendations In the use of the nonisothermal model to predict web tension, it is recommended that the temperature functions of the properties be determined by experimental measurement, especially, that ofthe Young's modulus. Caution should also be exercised in selecting the number of grid points to ensure precise prediction. Future investigation should be aimed at reducing the limitations of the model posed by the following assumptions. Assumption 2: The web material behaves purely elastically. Assumption 3: The web system is operating at steady state. Assumption 10: Viscous dissipation is insignificant compared to surfaced convection. The VEM model developed by Guan (1995) is capable of solving isothermal, mUltispan, viscoelastic problems under both steady state and transient conditions. The merge of the VEM model and the nonisothermal model will eliminate Assumptions 2 and 3. The resulting model can have the combined benefit to solve web handling problems with both thermal and viscoelastic effects. However, there are several changes that have to be made in the model development phase. 1, The significance of viscous dissipation (Assumption 10) needs to be studied. 2, A time dependent term should be included in the energy balance for unsteady 88  state, leading to a partial differential equation with two spatial variables and one time variable. 3. The force balance in the VEM model should be redeveloped to account for thennal strain. 4. Simplifications made under the assumption of isothermal conditions in the derivation of the VEM model should be reworked. 5. The temperature dependence of the material properties should be reflected in the VEM model. 89   REFERENCES 1. Aklonis, J. J., MacKnight, W. J., and Shen, M., 1972, Introduction to Polymer Viscoelasticity, WileyInterscience, New York. 2. Ames, W. Fo, 1965, Nonlinear Partial Differential Equations in Engineering, Academic Press, New York. 3. Back, E. L., and Andersson, L. 1., 1993, "The Effect of Temperature on Wet Web Strength Properties," Tappi, Vol. 76, No.5, pp. 164172. 4. Baird, D. Go, and Collias, D. 1., )995, Polymer Processing Principles and Design, ButterworthHeinemann Series in Chemical Engineering, Boston. 5. Balasubramaniurn, S., 1997, "Development of Interactive Software to Study the Effects of Viscoelasticity in Web Lines," Master's Thesis, Oklahoma State University. 6. Bicerano, 1., 1993, Prediction of Polymer Properties, Marcel Dekker, Inc., New York. 7. Bird, Ro 8., Stewart, W. E., and Lightfoot, E. N., 1960, Transport Phenomena, John Wiley & Sons, New York. 8. Bird, R. B., Annstrong, R. c., and Hassager, 0.,1987, Dynamics of Polymeric Liquids, Volume 1, Fluid Mechanics, Second Edition, John Wiley & Sons, New York. 9. Boresi, A. P., 1965, Elasticity in Engineering Mechanics, PrenticeHall, Inc., Englewood Cliffs, New Jersey 10. Boyd, 1. Po, 1989, Chebyshev and Fourier Spectral Methods, SpringerVerlag, New York. 11. Brandenburg, G., 1976, "New Mathematical Models for Web Tension and Register Error," Proc. 1, 3rd International IFAC Conference on the Instrumentation and Automation in the Paper, Rubber and Plastics Industries, Brussels, May 2426. 12. Brezinski, J. P., 1956, "The Creep Properties of Paper," Tappi, Vol. 39, No.2, pp. 116128.  90 1. _   13. Burnett, D. S., 1987, Finite Element Analysis: From Concepts to Applications, AddisonWesley, Reading, Massachusetts. 14. Byrd, Y. L., 1972, "Effect of Relative Humidity Changes During Creep on Handsheet Paper Properties," Tappi, Vol. 55, No.2, pp. 247252. 15. Campbell, D. P., 1958, Process Dynamics, John Wiley & Sons, New York. 16. Canuto, C., Hussaini, M. Y., Quarteroni, A., and Zang, T. A, 1988, Spectral Methods in Fluid Dynamics, SpringerVerlag, New York. 17. Carslaw, H. S., and Jaeger, J. C., 1959, Conduction of Heat in Solids, Second Edition, Oxford University Press, Oxford. 18. Chandrupatla, T. R., and Belegundu, A D., 1991, Introduction to Finite Elements in Engineering, Prentice Hall, Upper Saddle River, New Jersey. 19. Dugdale, D. S., 1968, Elements of Elasticity, Pergamon Press Ltd., Oxford 20. Dunn, H. S., 1969, "Equivalent Circuit Representation of Web Conveyance Systems," Journal of Applied Mechanics, June, pp. 316318. 21. Ferry, J. D., 1970, Viscoelastic Properties of Polymers, Second Edition, John Wiley & Sons, Inc., New York. 22. Finlayson, B. A., 1980, Nonlinear Analysis in Chemical Engineering, McGrawHill, New York. 23. Forsythe G. E., and Wasow, W.R., FiniteDifference Methods for Partial Differential Equations, John Wiley & Sons, New York. 24. Gao, 1., and Weiner, 1. H., 1992, "Computer Simulation of Viscoelasticity in Polymer Melts," Macromolecules, No. 25, pp. 13481356. 25. Gear, C. W., 1971, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, Englewood Cliff, New Jersey. 26. Gehlbach, L. S., Kedl, D. M., and Good, J. K., 1989, "Predicting Shear Wrinkles in Web Spans," Tappi Journal, August, pp. 129134. 27. Gottlieb, D. and Orszag, S. A., 1977, Numerical Analysis of Spectral Methods: Theory and applications, S.lA.M., Philadelphia. 28. Grenfell, K. P., 1963, "Tension Control on PaperMaking and Converting Machinery," Proc., IEEE Ninth Annual Conference on Electrical Engineering in 91  p:s the Pulp and Paper Industry, Boston Massachusetts, June 2021. 29. Guan, x., 1995, "Modeling of Viscoelastic Effects on Web Handling System Behavior," Doctoral Thesis, Oklahoma State University. 30. Harper, C. A, 1975, Handbook of Plastics and Elastomers, McGrawHill, Inc., New York. 31. Hartnet, J. P., 1992, "Viscoelastic fluids: A New Challenge in Heat Transfer," Transactions of the ASME, Vol. 114, May, pp. 296303. 32. Hauptmann, E. G., and Cutshall, K. A, 1977, "Dynamic Mechanical Properties of Wet Paper Webs," Tappi, Vol. 60, No. 10, pp. 106108. 33. Henderson, G. M., Hung, 1. C. and Googe, 1. M., 1979, "Automatic Tension Control for Winding Superconducting Coils," Proc. Southeast Conf. Reg. 3 Conf., Roanoke, Virginia. 34. Hilado, C. 1., 1969, Flammability Handbook for Plastics, Technomic Publications, Stamford, Connecticut. 35. Hindmarsh, A C., 1972, GEAR: Ordinary Differential Equation System Solver, Rep. UCID30001, Rev.2, Lawrence Livennore Lab., Livennore, California. 36. Hindmarsh, A C., 1973, GEARB: Solution of Ordinary Differential Equations Having Banded Jacobian, Rep. UCID30059, Lawrence Livermore Lab., Livennore, California. 37. Hindmarsh, A C. 1983, "ODEPACK, A Systematized Collections of ODE Solvers," In Scientific Computing, Stepleman, R. S., et al. (eds.), NorthHolland, Amsterdam, pp. 5564. 38. Hindmarsh, A c., 1987, "Brief Description ofODEPACK  A Systematized Collection ofODE Solvers," http://www.netlib.org/odepack. 39. Hopfinger, A J., Koehler, M. G., and Pearlstein, R. A, 1988, "Molecular Modeling of Polymers. IV. Estimation of Glass Transition Temperatures," Joumal of Polymer Science: Part B
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Modeling of Thermal Effects on a Single Span Web Handling System 
Date  19990501 
Author  Liu, Jian 
Document Type  
Full Text Type  Open Access 
Note  Thesis 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  MODELING OF THERMAL EFFECTS ON A SINGLE SPAN WEB HANDLING SYSTEM By JIAN LIU Bachelor of Science Oklahoma State University Stillwater, Oklahoma 1995 Submitted to the faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE May, 1999 MODELING OF THERMAL EFFECTS ON A SINGLE SPAN WEB HANDLING SYSTEM Thesis Approved: Thesis Adviser ii ACKNOWLEGEMENTS I wish to express my special appreciation to my thesis advisor Dr. D. Alan Tree for his invaluable guidance, inspiration and encouragement throughout both my undergraduate and graduate study at Oklahoma State University. Dr. Tree has been a true mentor who offers his most sincere advice in both my academic pursuit and life endeavor. I am also very grateful for the suggestions and support from the other members of the thesis committee, Dr. Martin S. High and Dr. Jan Wagner. Special thanks are due to the School of Chemical Engineering for furnishing me with the best knowledge and understanding of Chemical Engineering that the world can offer. Generous financial support in the past three years from the School of Chemical Engineering and the Oklahoma Center for Integrated Design and Manufacture is gratefully acknowledged. Finally, I would like to dedicate this thesis to my family. I am compelled to do so by their tremendous love, support and sacrifices. It is by God's grace that I was born in such a supportive and loving family. iii Chapter TABLE OF CONTENTS Page I. INTRODUCTION 1 II. BACKGROUND 4 2.1 Web Handling Models 4 2.2 Viscoelastic Effects in Web Handling 8 2.3 Temperature Dependence of Web Properties 11 2.4 Numerical Methods 19 2.5 Summary 21 III. MODEL DEVELOPEMNT 23 3.1 System Analysis 23 3.2 Assumptions 26 3.3 Governing Equations 28 3.3.1 Mass Balance 28 3.3.2 Force Balance 29 3.3.3 Energy Balance 32 3.3.4 Initial and Boundary Conditions 37 3.4 Nondimensionalization 38 3.5 An Analytical Solution for Constant Properties 40 3.6 Summary 45 IV. NUMERICAL METHODS AND TECHNIQUES 46 4.1 Selection of Numerical Techniques 46 4.1.1 Method of Lines and PDEONE 48 4.1.2 LSODE 51 4.1.3 Overall Algorithm Structure and Subroutine Arguments 53 4.2 Tension Calculation 59 iv  V. RESULTS AND DISCUSSION 62 5.1 Verification 63 5.2 Property Functions and Parameters 66 5.3 Results 70 5.3.1 Temperature Profile 71 5.3.2 Strains and Tensions 74 5.3.3 Grid Size 78 5.4 Parametric Studies 82 VI. CONCLUSIONS AND RECOMMENDATIONS 86 6.1 Conclusions 86 6.2 Recommendations 88 REFERENCES 90 APPENDIXES 97 APPENDIX AProgram Listings of the Main Program and User Defined Subroutines 98 APPENDIX BPhysical Properties of Several Web Materials 106 APPENDIX CParametric Studies of the Temperature Profile in the Thickness Direction 110 APPENDIX DEstimation of Viscous Dissipation 113 v LIST OF TABLES Table Page 1. Argwnents of the Three User Defined Subroutines 56 2. Arguments of LSODE 58 3. Test Case System Parameters and Dimensionless Groups 63 4. Fast Convergence of the Analytical Solution 64 5. Operating Parameters of a Typical System 68 D.I. Estimated Values of the Variables in the Energy Equation 114 vi UST OF FIGURES Rg~ P~ 1. Temperature Dependency of Polypropylene Properties 24 2. Thennal Strain Illustration 26 3. A Control Volume Cut from a Web Line (not to scale) 28 4. Elongation due to Defonnation and Thermal Expansion 30 5. A Web Span with Initial and Boundary Conditions 37 6. Node Point Layout of the MOL Example Problem 49 7. Layout of the Numerical Approach 54 8. Overall Algorithm Structure 55 9. Dimensionless Temperature Profile of Test Cases 65 10. Web System Dimension Summary 68 11. Isothermal Tension at Different Temperature 70 12. Dimensionless Temperature Profile of a Single Polypropylene Web Span 71 13. Temperature Profile ofa Single Polypropylene Web Span 72 14. Dimensionless Temperature Profile ofa PP Web Span over Thickness 73 15. Strains over a Range of Initial Temperatures 74 16. Tensions over a Range ofInitial Temperatures 76 17. Tensions over a Range of Ambient Temperatures 77 vii 18. The Nonisothennal Tension Prediction Error of Different Grid Sizes 80 19. The Thennal Expansion Prediction Error of Different Grid Sizes 81 20. The Effect ofPechlet Number on the Dimensionless Temperature 83 21. The Effect ofNusselt Number on the Dimensionless Temperature 84 c.l. Dimensionless Temperature Profiles at p=0.1 for Varioas Pe's III C.2. Dimensionless Temperature Profiles at p= 0.5 for Various Pe's III C.3. Dimensionless Temperature Profiles at P= I for Various Pe's 112 CA. Dimensionless Temperature Profiles at p= 1 for Various Nu's 112 vijj A b C l p D E EI F G h k L m NOMENCLATURE web crosssectional area one half ofweb thickness heat capacity at constant pressure, at constant volume constant pressure heat capacity of liquid state polymer diffusion coefficient Young's modulus elongational strain in transport direction tension shear modulus shear modulus in an uncrosslinked polymer in the rubbery plateau surface heat transfer coefficient thennal conductivity constant, mass flow rate endtoend length of a single repeat unit length of web power law coefficient ix M n N NBBrot NSGrot NPDE NEQ Nu p Pe q R Rc t T molecular weight of a repeat unit entanglement molecular weight (a) power law exponent, (b) number of grid points in transport direction number of vertices in hydrogensuppressed graph of repeat unit number of rotational degrees of freedom in the backbone of a repeat unit correction tenn for standard molar volume correlation total number of rotational degrees of freedom in the side group of a repeat unit Tg correlation parameter of twelve structure parameters number ofpartial differential equation(s) number of ordinary differential equation(s) Nusselt number pressure Pechlet number heat flux ideal gas law constant ratio of half web thickness over web length time temperature x u V, VI, V2 v w Xu ambient temperature temperature at web boundary temperature of the ith element (a) initial temperature, (b) reference state temperature glass transition temperature internal energy web velocity, first roller velocity, second roller velocity volume van der Waals volume web width thirteenth structure parameter for Tg Greek Symbols a (a) incline angle of web surfaces, (b) dimensionless thickness coefficient of volumetric thermal expansion (a) incline angle of web sides, (b) dimensionless length coefficient of linear thermal expansion zerothorder simple connectivity index zerothorder valence connectivity index firstorder simple connectivity index xi L1L E: y y or y = =(1) y 110 8 p a "C =(1) u Captions 1\ firstorder valence connectivity index solubility parameter difference operator length between two grid points in the transport direction total strain elastic strain thennal strain shear strain tensor shear strain rate tensor or convected time derivative of shear strain tensor magnitude (second invariant) ofthe shear strain rate tensor zeroshearrate viscosity time constant in the Maxwell model dimensionless temperature density nonnal stress viscous stress tensor contravariant convected time derivative of stress tensor Poisson's ratio vector indicator molar property xii ...... Subscripts o x, y, z at reference condition unless otherwise specified indicating X, y, zdirections, respectively Mathematical Symbols  'V I D gradient operator summation operator substantial derivative operator xiii CHAPTER I INTRODUCTION In the material manufacturing and processing industry, highly flexible, continuous thin sheets and strips are referred to as webs. Polymer films, paper, textiles and metal foils are all webs. Compared to large, discrete sheets of material, webs are easier to process, transport and less expensive to manufacture. Furthermore, web handling systems can be extensively automated to yield better product quality. Web handling is a widespread and essential process in the material manufacturing and processing industry. Because of the high speed of web lines, web handling poses many intriguing practical questions. However little systematic study was done before the inception of Web Handling Research Center (WHRC) at Oklahoma State University. The mission of the WHRC is to advance the knowledge base in technologies applicable to the transport and control of continuous strip materials through fundamental generic studies. Due to the increasing speed of webtransport systems, accurate tension control has been one of the major research areas at WHRC. Poor tension control may lead to wrinkling of the web, breakage or slackness, thereby limiting utilization of plant equipment, product quality, and profitability. A computerbased program, known as Web Transport System (WTS) (Shin 1991), has been developed for use in the analysis and design of openloop and closedloop control in webtransport systems. The initial 1 versions of WTS made some iC:~aUzdhai'~llmptionsfor the ~Y1Cbtmakri:lt,s~cS:lli:basjpure elastic behavior, no phase changes, no heat transfer, land no int.:rfaeiaima "IS,' tnmsl~r.J" The nonideal effec~s iIt the dmsron icontrobo.:f web lines inclJUde slj'1paf:.fb}1we~tra.M(.·"'·'Pl and a roller, web crosssectiond aIlIea change, composition chan~..;~.;rt>...:IJ.l1~',r~t'fJfechanr~ ~ and viscoelastic effects (Shin, 1991). The WHRC.has 'made a continuing effort to model the nGJ'ni dl iF efLet:9".in :v•.'.cb~~, transport systems and to incorporate these new models .u thl~ lat·~·stl,~':·rsi(','ll')llJf. WT~T;1S ;.. the first step to model nonideal effects, Guan (1995) developed a viscoelastic model (VEM) to study and predict tension due to viscoelastic behaviors. in WI6hhandling processes. In order thtisolate thte viscoelastic behavior and simplify'the pro~Lem",Guan continued to use a few of the ideal assumptioFls in Shin's work, including uniaxial tC'Iilsion:: and isotropic web materials. The results of Guan's work demonstn:.wd,the irTlplDTJtJodeJ,oL .'. viscoelasticity in w.eb:handling, as wellas,pointing out th .... ntfCossitt) fa~':oQd;yio1'hJ,:t~:r nonideal effects. This study was motivated by the need ·for it, nonisoth"rma,l moch:.:t,H't.J .:~~ address the issue of heat transfer. Heat transfer between the web and the environment is encountered in treatment and manufacturing processes such as cooling, drying, printing, and coating. For example, a frost line can be easily located in cooling sections where high temperature polymeric web materials undergo a phase change by losing heat to the environment. Since all material properties are functions of temperature, it is important to understand the thermaJo:n.:·,l effects on a web handling system in order to have accurate tension control. This thesis will concentrate on dealing with heat transfer in webtransport systems. The remainder ofthis thesis is divided in to five chapters. Chapter 11 includes tt,,,, ,!'; 2 relevant background and a literature review. Chapter III contains the development of the nonisothermal model from the constitutive equations. The systems of equations are then nondimensionalized and reduced for the purpose of general analysis. Nondimensionalized initial and boundary conditions are also presented for the governing equations. Chapter IV describes the numerical methods used to solve the resulting nonlinear partial differential equation, and develops the tension formulation. Chapter V presents the numerical simulation results and discussion. Finally, Chapter VI concludes the thesis, and comments on the direction of future studies. 3 CHAPTER II BACKGROUND Web handling is a significant sector in the modem materials processing industry, and has attracted some research interest. However, little has been done on modeling the impact of nonideal effects on tension control due to the complexity of developing rigorous models. This chapter first reviews the development ofweb handling models, which is followed by a discussion of viscoelastic effects in web handling. Section 2.3 summarizes the available methods to determine the temperature dependency of web properties. Finally, a brief discussion on applicable numerical methods for partial differential equations is given in Section 2.4. 2.1 Web Handling Models The earliest systematic study of web handling systems can be traced to the late 1950s. Campbell (1958) examined the idealized web propulsion, tension control and guidance of web materials in a free span. Campbell proposed a tension model based on the assumption that the web material obeyed Hooke's law. In an attempt to capture viscoelastic effects, the free span tension formulation was then adjusted to include an empirical internal coefficient of viscous friction. This rather simple adjustment did not 4 satisfactorily account for viscoelasticity, and no other attempt was made to capture nonideal effects. In addition, Campbell did not account for the well known "tension transfer" to the span of interest, which results from the material being in tension in preViOUS span. Assuming that web materials obeyed Hooke's law, Grenfell (1963) incorporated the tension transfer effect into his formulation and modeled the interaction of adjacent spans in multispan systems. Grenfell's model was capable ofdescribing dynamic conditions such as a step change and a sinusoidal variation in the speed of the end roller. Grenfell found that web tensions did not instantaneously respond to a change in speed. A time constant was used to measure the delay in tension response. Furthermore, Grenfell pointed out that slippage and change of moisture content in paper sheets could upset the web tension. Grenfell's work proposed feed back control strategies to solve these problems. Brandenburg (1976) reviewed tension studies of web lines, and developed a mathematical model to answer the increasing need for processing automation. Brandenburg continued to assume that the web behaved elastically and that no mass transfer occurred in and out ofthe web. The Brandenburg model successfully accounted for the adhesion and slippage of the web in the contact region with a roller due to different velocities in adjacent spans. Brandenburg's isothermal study showed that the tensile force, stress and strain were position independent in a free span. This conclusion is incorrect in processes affected by nonideal conditions such as viscoelasticity, heat transfer, and mass transfer. Furthermore, Brandenburg modeled the effect of changes in temperature both in 5  the span and in the contact region. Brandenburg introduced the thermal strain to the elastic strain to obtain the total strain. However, the temperature dependency of web properties was not included in this model. Heat transfer between the web and the environment was modeled using a thermal time constant and an average temperature at a given position. Brandenburg's model did not solve for the temperature profile within the web, which is needed in determining local properties in a rigorous mathematical model. Taguchi et al. (1985) realized that many processing problems, e.g. web breaks, wrinkles, cutoff errors and unstable runs, could be solved by adjustment of web tension. Taguchi's experimental work also showed that the tension transfer was significantly affected by the slippage between the web and the driven roller. A computer simulation, which assumed elastic behavior of the web, was developed to study the optimum design and control system for web tension. Shin (1991) developed a model which included seven primitive elements, i.e. unwinding roll, free span, driven roller, idle roller, dancer roller, nip roller and winding roll. Combinations of these seven elements can be used to describe a complex web handling operation. This idea resembles the "unit operation" concept in chemical engineering processes. Shin reviewed and tried to summarized the impact of nonideal conditions of web material, such as temperature change, moisture change and web viscoelasticity. Shin's unified model for a multispan system included the effect of slippage, temperature variation and moisture variation. Viscoelastic behavior was not considered. Shin argued that closed loop control is necessary to solve the tension transfer problem due to slippage. In a closed loop scheme, the web handling systems are treated as interconnected 6 subsystems made up of a roller and a web span. All subsystems speeds are measured and coordinated simultaneously by a master controller applying the unified model. Although complicated, the close loop control concept has been widely discussed in the field of web handling. The major drawback of the closed loop control is the difficulty in obtaining accurate online tension measurements. Shin's approach to include the effect of temperature variation ofthe web required that the web temperature be known and input to the model. This approach did not capture the essence of heat transfer within the web and between the environment and the web. Following Shin's work, Reid and Lin (1993) studied the dynamic behavior of multispan systems. They utilized simple fixedgain PID controllers in their simulation examples, and concluded that variablegain controllers were needed to produce good dynamic behavior over a broad range of startup conditions. There are other models developed for different aspects of web handling systems. Although the importance of nonideal effects has been known for decades, many models neglected them due to the complexity of introducing even one of these nonideal effects. In addition to the models already discussed in this section, there are similar elastic models in studies of multispan systems (Dunn, 1969; Soong and Li 1979; Young et a1., ]989a; Young et a1., 1989b; Parant et a1., 1992), tension control (Martin, 1973; Henderson et aI., 1979), tension measurements (Horst and Negin, ]992), tension oscillations (Veits et aI., 1981), lateral dynamics (Shelton and Reid, 1971 a; Shelton and Reid, 1971 b), winding mechanisms (Pfeiffer, 1966; Pfeiffer, 1968; Rand and Eriksson, ]973; Pfeiffer, 1977), and wrinkles (Gehbach et a1., 1989). A few viscoelastic studies can also be found in winding mechanics (Tramposch, ]967; Mukherjee, 1974; Lin and Westmann, 1989; 7  Kalker, 1991), and single span behavior of paper handling (Hauptmann and Cutshall 1977). Among these models, Tramposch (1967) and Parant (1992) incorporated temperature effects in their models, but made no attempt to find the temperature distribution within the web. Because ofthe highly temperature sensitive viscoelastic properties, Guan (1995) pointed out the importance of obtaining the temperature distribution to accurate tension modeling and control. There are also several papers dealing with the moisture and temperature effects on web materials. However, these studies were not in the area of modeling web handling systems. Morland and Lee (1960) worked on the stress analysis for linear viscoelastic materials with temperature variation. Byrd (1972) and Roisum (1992) both investigated the moisture effect on webs. LawaI (1989), Losi and Knauss (1992) and Back and Andersson (1993) studied the temperature effects on weblike materials. 2.2 Viscoelastic Effects in Web Handling When a viscoelastic material is deformed, it exhibits both elastic and viscous components. The elastic component is manifest by the recoverable strain energy. The viscous component results in an irrecoverable energy dissipation. The ratio of viscous to elastic response of a material depends on its microstructures and the processing conditions, which include temperature, pressure, geometry of deformation, type of force, prehistory and the time interval of deformation. Most web materials exhibit viscoelastic behavior. Many commercial web handling systems are designed to operate in the region where the materials exhibit 8 primarily elastic behavior to avoid complications in the tension control. However, there are many other processes or sections of processes in which the correct modeling and control in viscoelastic web materials are essential to product quality and production management. For example, the newsprint systems are heavily dependent on the viscoelastic properties of the wet paper webs under dynamic conditions. The effect of moisture is very significant to the wet paper viscoelasticity due to the weak hydrogen bonds between the crossing paper fibers (Brezinski, 1956). The environmental temperature detennines the maximum amount of moisture the air can hold, and thus, affects the relative humidity of the processing environment. The relative humidity of the environment, in turn, affects how much and how fast the moisture of the wet paper evaporates. Therefore, three nonideal effects, viscoelasticity, temperature and moisture are encountered simultaneously in the wet paper web handling systems. For other polymeric materials, the viscoelasticity will generally increase as the processing temperature increases (Ferry, 1970; Aklonis et aI., 1972; Ward, 1983). Viscoelastic models have been under development for more than a century. There are currently three types of rheological equations of state to describe the viscoelastic behavior, i.e. the Generalized Newtonian Fluids (GNF), the linear viscoelastic models and the nonlinear viscoelastic models. However, no single rheological equation of state can universally predict the viscoelastic responses of all known viscoelastic materials in various applied situations (Spriggs et aI., 1966; Tanner, 1983; Bird et aI., 1987). The GNFs are the simplest constitutive equations. The stress and strain tensors are related by a viscosity function, which is commonly chosen to be a power law relationship to account for nonNewtonian viscosity. The GNFs are useful in modeling 9 behavior under steady shear flows conditions.. The GNFs' inability to predict elastic effects limits their application in web handling systems where the elastic component is usually nontrivial. The General Linear Viscoelastic Fluids are a large group of rheological equations of state, which linearly combine Hooke's law and the Newton's law of viscosity. The idea of this type of models is to simulate the viscoelastic response using various parallel and series arrangements of a spring and a dashpot. The Maxwell model is a typical example that uses a series combination, and results in the following expression: '\ a . "( + Al "( =TJo y. = at = = (2.1 ) Th.e Voigt model, the generalized Maxwell model and the Jeffreys model are some of the other common models in the linear viscoelastic group. Guan (1995) pointed out that the linear viscoelastic models lack material objectivity when applied to systems in motion. These models predict different results for the same moving system depending on the coordinates chosen for the model analysis. The problem with material objectivity can be traced to the partial derivatives of the linear viscoelastic models. Since all web handling systems involve translation motion, Guan did not select a linear viscoelastic model. The nonlinear viscoelastic models overcome the material objectivity problem of the linear models by introducing convected derivatives to replace the partial derivatives. Guan (1995) selected the WhiteMetzner model (Bird et aI., 1987) for the study of viscoelastic effects on web handling systems: "( + TJ(Y) '[ = TJ(Y)y , = G =(1) =(1) 10 (2.2) D where (2.3) The subscript (1) denotes the upper convected derivative. The WhiteMetzner model is a nonlinear Maxwell model with a modification to include the nonNewtonian viscosity function. The WhiteMetzner model allows for stress relaxation through viscous deformation; preserves material objectivity, independence of results on coordinate system; provides the flexibility to model diverse materials such as polymer and paper; and has sufficient mathematical simplicity to allow for a solution (Guan, 1995). 2.3 Temperature Dependence ofWeb Properties The physical properties of web materials are functions of temperature. There are three ways to find the propertytemperature relations. The first, and most accurate way, is to experimentally measure the property functions of a particular material of interest. However, measuring all the needed property functions is not efficient in any preliminary design stage. The second way is to find published data on similar materials. Tremendous amounts of experiments have been done to investigate and document the physical properties of polymeric materials. Mark (1996) has organized the literature information on the physical properties of polymers. For example, the density of atactic poiypropylene as a function of temperature for various polymers is given along with the applicable temperature range and references: (2.4)  The temperature range for this correlation is 80ae  120D e. Eq. (2.4) was regressed from 11 (2.5)  data published by Rogers (1993). Although not given over a range of temperature, the thermal conductivity and heat capacity of polypropylene can also be found under other sections of Mark's work. Other properties such as the elastic modulus and the power law coefficients are not available through Mark's work. Heat capacity and thermal conductivity data can be found in Steere's work (1966) among many other published works (Harper, 1975; Hilado, 1969). The elastic modulus can be obtained from Reding (1958). This raises the problem of searching through a large amount of literature for the data in the desired form. Since the discovery of polymeric materials, many models have been developed to predict polYmer properties. As polymers with relatively simple repeat unit structure have reached their limits, more new polYmers are being synthesized for different applications. Some models are capable of making predictions even before a polYmer is synthesized. One comprehensive and very successful prediction is the group contribution techniques in the qualitative structureproperty relationships (QSPR). Van Krevelen (1976, 1990) published the classic textbook that contains a compendium of useful QSPR relationships for polymers. The group contribution techniques conceptually break down a molecule into small fragments or groups of atoms. The total value of a molar property is the summation of all the group contributions to the property of interest. The idea is illustrated by the following equation. (Total Property) :::::: I("additive" group contribution) + I("constitutive" structural terms). In Eq.(2.5), the second summation over structural terms only has to be included for certain properties where the assumption of strict addition of group contributions results in 12 p  large errors because the environment of a given group also plays an important role in the contribution of the group to that particular property. The values ofindividual group contributions can be estimated by fitting the observed values of the properties of interest from experimental data on other polymers containing the same types of fragments. Therefore, the group contribution techniques are essentially empirical and rely on statistical analysis of data. This explains the successful prediction of some properties of relative simple structure to within 1% of experimental values. In spite of its usefulness, the group contribution techniques have a few inherent limitations. The most apparent and important limitation is the need for a substantial body of experimental data for a given property in order to derive reliable values for the contributions of a large number of molecular fragments. Attempts were made to bypass the need for experimental data with theoretical models of molecular interactions. However, this approach requires extensive amounts of timeconsuming force field and quantum mechanical calculations (Hopfinger et al., 1988). Seitz (1993) developed a set of correlations to simplify the group contribution techniques. This correlation allows the prediction of a large number of "derived" physical properties of polymers from just five "fundamental" properties: molecular weight of repeat unit, endtoend length ofrepeat unit in full extension, Van der Waals volume of the repeat unit, cohesive energy, and a parameter related to the number of rotational degrees of freedom of the backbone of a polymer chain. This approach gives reasonably accurate predictions, provided the five fundamental properties can all be estimated accurately. The major drawback of Seitz's simplification is the inability to estimate the cohesive energy (Bicerano, 1993). 13 •  Bicerano (1993) suggested an alternative approach, called the topological technique. The topology is the pattern of interconnections between atoms. Since only about a dozen nonmetallic atoms (e.g. C, Hand 0) made up a large number of common macromolecules and there are only certain spatial bonding arrangements of these atoms. The development and fonnalism of "connectivity indices" can be used to mathematically predict the properties ofmacromolecules. There are significant improvements, especially in model simplicity, over other contemporary methodologies. The detailed treatment and application of the connectivity indices to different physical properties are beyond the scope of this thesis. One can find the relevant derivations in the book published by Bicerano (1993). Bicerano first obtained the four connectivity indices from experimental data. These four connectivity indices are the zerothorder simple connectivity index ox., the zerothorder valence connectivity index 0x.v , the first order simple connectivity index lX, and the first order valence connectivity index IXV. The values of these connectivity indices are tabulated for many common polymers by Bicerano. For example, the connectivity indices for polyethylene are 1.4142, 1.4142, 1.0000 and 1.000 respectively, and for polypropylene 2.2845,2.2845, 1.3938 and 1.3938. The standard properties at room temperature are then expressed as functions ofthe connectivity indices or other standard properties. The properties at any temperature are calculated from a correlation based on the standard properties. The coefficients of these correlations are obtained empirically, i.e. by curve fitting experimental data. The results of the topological approach to find the temperature correlation of the needed properties are presented below, followed by a brief discussion of a few limitations of the topology 14 technique. The temperature dependence of the molar volume of polymers are given in two sets of correlations. For polymers with Tg ~ 298 K, 1.42 Tg + 0.15 T V(T) =V(298 K) 1.42 T g + 44.7 ' (for T s Tg and Tg ~ 298 K); (2.6) 1.57 Tg +0.30 (T  Tg ) V(T) = V(298 K) , (for T > Tg~ 298 K). (2.7) 1.42 Tg + 44.7 For polymers with Tg s 298 K, (for T < Tg < 298 K); (2.8) Vet) = V(298 K) [1 + a r (298 K) (T  298)], (for T ~ Tg and Tg < 298 K), (2.9) where aT (298 K) = ( . ) , 298 + 4.23 Tg (for Tg < 298 K). (2.10) The molar volume at room temperature V(298 K) can be determined by a correlation of the four connectivity indices and a correction term NMy. V(298 K) =3.642770 Ox +9.79897 °X v 8.542819 IX + 31.693912 IXv + 0.978655 N MY' The value ofNMv is 0 for both polyethylene and polypropylene. (2.11 )  The glass transition temperature of different polymers are fitted into a correlation of four parameters, which is tabulated in Table 6.2 of Bicerano (1993). 15  N'Ig Tg =351.00+5.63 8+31.68N 23.94 x 13 . (2.12) The four parameters are the solubility parameter 8, the correlation parameter of twelve structure parameters NTg, the number of vertices in the hydrogensuppressed graph of the repeat unit N, and the thirteenth structure parameter Xl3. The tabulated values of parameters 8, NTg, N and XI3 are 15.9, 16, 2, 0.0, respectively, for polyethylene; and 17.0, 20,3, 0.0, respectively, for polypropylene. The heat capacity of solid polymers as a function oftemperature is commonly extrapolated based on the heat capacity of the polymer at a reference temperature. Bicerano used the heat capacity at room temperature, which can be predicted using a correlation of four parameters. Only the correlations for "liquid" state polymers were chosen because most polymeric materials are in the rubbery state in web handling applications. C~ (T) = C~ (298 K) [1 + 1.3 x 103 (T  298)]. C~ (298 K) =8.162061 Ox + 23.215188 Ox v + 8.477370 NBBrOl + 5.350331 N SGroI' (2.13) (2.14) where NBSro! is the number of rotational degrees of freedom in the backbone of a repeat unit, and NSGrot is the total number of rotational degrees of freedom in the side groups of a repeat unit. For polyethylene and polypropylene, the tabulated values ofNBBro! arld NSGrol are 2.0, 0 and 2.0, 1. Other than the above mentioned characteristic material properties (e.g. Tg, V, Cp), all mechanical properties of polymers (e.g. the Young's modulus) depend on the 16 > specimen properties which include crystallinity, orientation, crosslinking, thermal history and the type of force applied. This means predictions of mechanical properties are very specimen dependent. Experimental measurement is recommended to obtain reliable correlations of mechanical properties as a function of temperature. However, Bicerano's predictions are presented as a "quick and dirty" estimation for the Young's modulus E of uncrosslinked polymers under small strain in the temperature range ofthe rubbery plateau. The shear modulus ONo(T) of an uncrosslinked polymer in the rubbery plateau, is usually assumed to be determined by the physical interactions caused by entanglement between polymer chains. With the introduction of the entanglement molecular weight Me, GNO(T) can be expressed as a function of the gas constant, density, temperature and Me: G~(T) = p(T) R T. Me Me can be very roughly estimated as: N V M Me =1039.7 + 1.36411 X to23 BBrot 3 w 1m (2.15) (2.16) where NBBrot is the number of rotational degrees of freedom in the backbone of a repeat unit; M is the molecular weight of a repeat unit; Vw is the van der Waals volume; and 1111 is the endtoend length of a single repeat unit. The estimated values of these parameters are 2, 28 g/mole, 20.5 cc/mole, 2.52 x to8 em, respectively, for polyethylene and 2,44 g/mole, 30.75 cc/mole, 2.52 x to8 em, respectively, for polypropylene. For polymers in the rubbery region (T> Tg + 300 e), the Poisson's ratio, v, is very close to 0.5. The Young's modulus relationship with the shear modulus and the 17 • Poisson's ratio results in the following simple expression for uncrosslinked polymers in the rubbery region. (2.17) The thermal conductivity of amorphous polymers can be correlated from the value at Tg. (T> Tg). (2.18) lfthe thermal conductivity at any temperature is known, k(Tg) can be calculated using Eg. (2.17) and, hence, at all other temperatures. Therefore, k(298 K) is used as the reference value, and further estimation functions are needed. II [C~(298 K)] [ UR ]3 k(298 K) =5 ·10 " V(298 K) . V(298 K) , (Tg <298 K). (2.19) The molar volume and standard liquid heat capacity that appear in Eq. (2.19) are given by Eq. (2.6) to Eq. (2.9) and Eg. (2.13). The new parameter UR is called the molar Rau function and has units of cmIO /3/(sJl3mole). The estimated value of UR is 1760 for polyethylene and 2730 for polypropylene. Bicerano (1993) did not give correlations for the power law coefficients of the nonNewtonian viscosity, i.e. m and n in Eq. (2.3). Baird & CoUias (1995) have suggested the following empirical relationship between temperature and the power law coefficients. (2.20) (2.21 )  18  The reference parameters mOand nO, and the constants Band C have to be determined from experimental data. There are two limitations of the Bicerano correlations for web materials. These correlations are mainly developed for isotropic, amorphous, linear polymers. The effects of polymer crystallinity, orientation, long chain branching and crosslinking on polymer properties are not incorporated. Secondly, the connectivity indices and the correlations are obtained empirically by considering the correlation of the topological features of molecules with the properties of interest. However, the topological approach does provide some improvement over the group contribution techniques (Bicerano 1993). Polymer properties are very structure dependent. Collecting experimental data on the specific polymer of interest and curve fitting the experimental data is always the most accurate way to estimate the temperature dependence of the properties. Published data of similar specimens will provide insight and a relatively good estimate. If both of these methods are not plausible or a quick estimate is needed, caution should be exercised when using models for prediction; because no model can account for all the structural details of a particular specimen without losing generality. 2.4 Numerical Methods The study of a nonisotherrnal web transport system requires invoking the energy balance that usually consists two or more independent variables (position and time). The conduction tenus in the energy balance are second order partial derivatives. Since the properties are temperature dependent, the governing equation and the boundary 19  conditions are often highly coupled. Analytical solutions to this type of partial differential equation can usually be accomplished only with over simplification or rough approximation in which the accuracy is often lost. Considering the efficiency of modem age computing, many researchers have chosen the numerical approach. Many numerical techniques have been successfully developed to solve a wide variety of partial differential equations. Most techniques fall into the category of the finite difference methods. There are also other methods such as finite element, Monte Carlo, spectral, variational, series expansion and eigenfunction expansion methods. Each of these methods has its strength in a certain area of applications. For example, the finite element methods are popular in the field of solid mechanics and structural engineering due to its considerable freedom in setting the location of computational elements in highly irregular geometry (Strang and Fix, 1973; Burnett, 1987; Chandrupatla and Belegundu, 199 I). The spectral methods converge more rapidly than finite element methods, but only work well for very regular geometry with smooth functions (Gottlieb and Orszag, 1977; Canuto et aI., 1988; Boyd, 1989). Many specific techniques also work well with a particular type of physical problem. Patankar (1991) has developed a computer program, CONDUCT, for the analyses of heat conduction and duct flow heat transfer. In spite of these developments, the numerical solution of physically realistic, nonlinear partial differential equations is a complicated and highly problem dependent process which usually requires a researcher to undertake the difficult and time consuming task of developing a computer program. In order to eliminate the expensive and time consuming effort to solve nonlinear 20 > partial differential equations, Madsen and Sincovec (1975a, 1975b) presented a powerful and generalized software interface, called PDEONE. Based on the method of lines, PDEONE discretizes the space dimension and converts onedimensional partial different equations into a system oftime dependent ordinary differential equations, which can be reliably solved by a robust integrator. Melgaard and Sincovec (1981a, 1981 b) later developed another software interface, PDETWO, which is dedicated to solve twodimensional partial differential equations in rectangular coordinates. Both PDEONE and PDETWO require a robust integrator for initial value problems of systems of ordinary differential equations. The authors of these interfaces suggest the ODE integrator, GEARB, developed by Hindmarsh (1972,1973). In this study, an upgraded version of GEARB, LSODE (Hindmarsh, 1983) was used as the integrator for the systems of ordinary differential equations obtained from semidiscretization by PDEONE. Section 4.1 explains why these numerical methods were chosen and how the numerical routines function with each other. 2.5 Summary At the beginning ofthis study, the nonideal effects on modeling web handling systems were beginning to be addressed. Guan (1995) has completed the modeling of viscoelastic effects on multispan web handling systems. However, there were no nonisothermal models to capture the temperature distribution within the web material under nonisothermal conditions. Existing nonisothermal models assumed an average temperature variation, and dealt with the thermal effects separately from the viscoelastic 21   effects. A rigorous nonisothermal model, therefore, is needed to supplement Guan's rigorous viscoelastic model. The desired model should include the temperature dependency of web physical properties, and show how the temperature distribution affects the web tension. There are also difficulties in solving the complex energy equation, which is a highly nonlinear partial differential equation. These situations establish the background for this study. 22 a   CHAPTER III MODEL DEVELOPMENT This chapter contains the model development for temperature effects in web lines. In Section 3.1, the thermal effects on a single span web handling system are analyzed. The related assumptions for the single span system are given in Section 3.2. In Section 3.3, a set of general governing equations and the boundary conditions are developed. Section 3.4 presents the derivation of the dimensionless form of the governing equation and boundary conditions. In Section 3.5, an analytical solution with constant material properties is derived for the verification of the numerical results and parametric studies of the nonisothermal model. 3.1 System Analysis Heating and cooling in web transport systems are characterized by heat transfer between the environment and the web line by surface convection. Heat transfer between the environment and the web line causes a temperature gradient in the web material. The temperature gradient affects the tension in the web line in two ways: (1) by affecting the physical properties and (2) by inducing a thermal strain. For example, the Young's modulus is very sensitive to temperature. Since the tension is calculated using the 23  Young's modulus, temperature variation can significantly impact tension. In addition, some moderately temperature dependent properties, such as heat capacity and thermal conductivity, can also affect the process of heat transfer. Figure 1 shows the temperature dependency of some physical properties of polypropylene. The Young's modulus (E) is so sensitive to temperature that E must be plotted in log scale to fit on the same graph with the other properties. At high temperature (> 120°C), E is about one order magnitude smaller than at room temperature. 8.0·., p x 102 (kg/m 3 ) 8.0 • • • • • • • • • I I • • • • I I I • • • • • I I I • .....,. 7.0 4.0 log(E(psi)) 3.0 80 100 120 140 160 Temperature 20 40 60 2.0 +~_ __+____+_____l o Figure 1. Temperature Dependency of Polypropylene Properties* * Density, Mark (1996); Young's modulus, Reding (1958); Heat capacity and thermal conductivity, Steere (1966). 24 l The tension level is expected to increase drastically as a polypropylene web cools. The other physical properties are the density (p), the constant pressure heat capacity (Cp), and the thennal conductivity (k). The increase in Cp with temperature implies that polypropylene is more resistant to thennal changes at high temperatures. The increase in k with temperature indicates that polypropylene has greater thermal conductivity at higher temperature due to the increasing freedom of movement of polymer chains. The density is the least sensitive property among the four properties. However, the effect of a small density change can be a very significant change in the web tension level through the generation of a thermal strain. A decreasing p with temperature means polypropylene shrinks when cooled (negative thermal strain), and expands when heated (positive thennal strain). Temperature change directly induces a thermal strain in a constrained web even in the absence of stress (Shin 1991). Figure 2 illustrates the idea of thermal strain. lJndeT no constraint, the material is free to shrink or expand. In the unconstrai ned case of Figure 2, the thermal strain is Et, but there is no tension. If constrained at both ends, the material contracts, and experiences an increase in tension when the temperature decreases. The total strain in the constrained material after the temperature change is the sum of the thermal strain and the initial mechanical strain before the temperature change. If the constrained material is heated, it expands and relaxes the tension level. If the thermal expansion is greater than the initial mechanical strain, being unable to support compression, the web material sacks. 25  Unconstrained Constrained I ~~~ : I i : I f1~. Z I High Temperature Low Temperature Figure 2. Thermal Strain Illustration For a highly viscous web material, heating due to viscous dissipation can be significant during the stretching. However, since most web systems have small draw ratios (i.e. small velocity difference between adjacent rollers), the dissipation effect is considered negligible (Guan, 1995). Therefore, the tension of a nonisothermal web transport system is produced by thermal expansion/contraction in the form of thermal strain, and mechanical drawing in the form of elastic strain. 3.2 Assumptions The thermal effects of a nonisothermal system can be modeled under a set of generic operating conditions. These operating conditions are presented as the following assumptions. 26    Assumptions: 1. The thickness of the web is much smaller than the width and length. 2. The web material behaves purely elastically. 3. The web system is operating at steady state. 4. The web is under uniaxial tension. 5. The incline angles are very small. 6. No shearing deformations exist in the planes perpendicular to the machine direction. 7. No slip between the web and the rollers. 8. The web is isotropic. 9. Inertial, gravitational and air traction forces are negligible. 10. Viscous dissipation is insignificant compared to surface convection. II. The heat transfer environments above and below the web are identical. This study further assumes that the web material is purely elastic due to the complexity in dealing with both the thermal and viscoelastic effects. Forces such as inertial, gravitational, air traction, and longitudinal tension that can potentially affect the web tension have also been neglected. Finally, this study focuses on single open span systems. Under steady state conditions, multiple span systems can be treated as a combination of single span systems due to the assumption of elastic behavior. 27 a  3.3 Governing Equations The governing equations for the single span were derived from the mass balance, force balance, and the energy balance. 3.3.1 Mass Balance Figure 3 shows that the control volume cut from a web line has a wedge shape with small inclined angles as stated in Assumption 5. Q w y x I1z .I I z (Top View) i t. (\ m?_~·z (Side View) Figure 3. A Control Volume Cut from a Web Line (not to scale) The equation of continuity of any control volume system states that the rate of mass accumulation is equal to the rate of mass going in minus the rate of mass going out. Hence, the mass balance of this control volume can be written as: 28 ....... • p(z, t)A(z, t)v z(z, t)  p(z + I!.z, t)A(z + I!.z, t)vz(z + I!.z, t) I!.z =[p(z, t + ~t)A(z, t + L1t) p(z, t)A(z, t)]. 6t (3.1 ) a A (= 2bw) is the crosssectional area of the web in the xy plane. Dividing the entire Eq. (3.1) by I!.z and taking the limit as I1z and 6t approach zero gives (3.2) Under steady state conditions (Assumption 3), all variables are independent of time. The right hand side ofEq. (3.2) vanishes. A simple final form ofthe mass balance results: pAvz =km =constan t. 3.3.2 Force Balance From Assumption 4, the web system is under uniaxial tension in the machine (3.3) direction represented by the zaxis. In other words, the stress components in both the x and y directions are zero. Only the stressstrain relationship in the z direction is considered. (3.4) Consequently, the stressstraintension relationship is straightforward. Fz = Acr z = AEE e . (3.5)  For the nonisothermal case, the strain must be reexamined. Figure 4a illustrates a material that is originally in an unstrained/unstressed condition. Figure 4b shows the 29   same material under an extensional force, and Figure 4c shows the material under an extensional force after being heated. Figure 4 illustrates that the elongation results from both the elastic defonnation due to the applied force and the thermal expansion. In the nonisothermal case, the strain will be based on the total elongation regardless of the cause. Figure 4. Elongation due to Elastic Deformation and Thermal Expansion In the case of an open span in a web handling system, Y 2 YJ t ::.......:....=E1.=E c +E. (3.6) The thermal strain for an isotropic material will be pure elongation or contraction with no shear component (Crandall, 1959). Therefore, the thennal strain can be written as: yttl 0 xy = Yyz = Yxz = , 30 (3.7) (3.8) where ~l is the coefficient of linear thennal expansion, £t is the thennal strain, l is the thennal shear strain, !1T is the change of temperature. For isotropic materials, the coefficient of linear thennal expansion equals onethird the coefficient of volumetric thennal expansion (at). (3.9) By definition, the coefficient of volumetric thennal expansion is the fractional rate of change of the specific volume as a function of temperature. Since specific volume is the inverse of density, the coefficient of volumetric thennal expansion coefficient can be related to density, which is a function of temperature. a' ~ M~) ~H:) , p p where V is the specific volume, p is the density, (3.10) and a subscript p is used to denote that the partial derivative is taken at constant pressure. Thus, given the density as a function of temperature for any web material, one can easily compute the coefficient of linear thennal expansion. The density of web materials can generally be fitted to a power series from experimental data bz peT) =A + BT +CT2+... , 31 (3.11 )  with T representing the temperature. By combining Eg. (3.9) through (3.11) and carrying out the partial differentiation, the coefficient of linear thermal expansion can be expressed in terms oftemperature: c I 1 B+2CT+... ~ = 2 3A+BT+CT +... (3.12) Expanding Eg. (3.6) in the machine direction using the stressstrain relationship, Eq. (3.4), and the linear thermal expansion coefficient (~t), Eg. (3.7), gives 0 E =_z+Rt~T z E I' . Combining Eq. (3.5) and Eg. (3.13) yields F E =_z +Rt~T. Z AE I' (3.13) (3.14) Since it is difficult to measure the crosssectional area of the web, Eq. (3.3) is used to substitute for A. F pv E =_Z__Z +~t~T Z k E m The web tension can be calculated if the relevant material properties and the web (3.15) b temperature profile are known. Eq. (3.15) will be formulated to a summation form later in Chapter IV due to the discrete temperature distribution of numerical solutions. 3.3.3 Energy Balance The energy balance of the control volume depicted in Section 3.3.1 is the most essential governing equation of this modeling study. The beginning is the equation of thermal energy presented by Bird, et al. (1960). 32   ps pDD =(v.q)p(v,v)(1::Vv). Ot (3.16) IQ The left hand side ofthe equation represents the rate of gain of internal energy per unit volume. The tenns on the right hand side of the equations are from the following energy contributions in their respective order: (l) the rate of internal energy input by conduction per unit volume) (2) the reversible rate of internal energy increase per unit volume by compression, and (3) the irreversible rate of internal energy increase per unit volume by viscous dissipation. The internal energy ofEq. (3.16) needs to be modified in terms of heat capacity and temperature. Bird et a1. (1960) suggested rewriting internal energy and specific volume by applying the chain rule. ...  (aD)  (aD) [ (ap) ]   dU = ay TdV + aT V dT =  P +T aT v dV + CydT )  (av) (aY) dV = aT pdT + ap T dp . Since most web transport systems operate under atmospheric pressure, web materials are considered under constant pressure. Eq. (3.] 8) becomes (~~) = ~~. p Combining the definition of heat capacity at constant pressure and Eq. (3.19) gives   (aY) (ap)  dY(ap) C =C +T   =C +T _. . p v aT or  v dT or  p v v Multiplying dT through Eq. (3.20) and rearranging the resulting equation yields 33  (3.] 7) n.18) (3.19) (3.20) .. (3.21 ) • Substituting Eq. (3.21) into Eq. (3.17) eliminates the tenns associated with heat capacity at constant volume. (3.22) Eq. (3.22) is rewritten in the fonn of substantial derivatives. Density is then multiplied on both sides of the equation. The continuity equation of any control volume can be expressed as: Dp _ =p(V·v). Dt Thus, the substantial derivative of specific volume can be transformed using the continuity equation (Bird, 1960). DV o(X) I Dp p=p ==(V·\i). Dt Dt p Dt Eq. (3.23) can be further simplified to DO  DT p=p(V·\i)+C p. Dt p Dt (3.23) (3.24) (3.25) (3.26) Finally, substituting Eq. (3.26) into the Eq. (3.16) yields a useful form of the equation of thennal energy. pCp DT =(V .q)  (1::V'\i). Dt (3.27)  This is the equation of thermal energy for a material at constant pressure. Eq. 34 (3.27) needs to be further simplified with the assumptions listed in Section 3.2. First, the substantial derivatives and the vector operations are expanded in a rectangular coordinate system.  (Of Of Of OfJ Pc +v +v +v  P at X ax Yay z8z (3.28) The viscous dissipation terms can be neglected due to small velocity variation in any directi<;>n (Assumption 10). All the terms associated with shear stress on the right hand side of the equation disappear (Assumption 6). Moreover, the velocities in x and y direction are negligible. Under steady state condition (Assumption 3), Eq. (3.28) is reduced to (3.29) The heat transfer area available for surface convection on both the top and bottom of the web, parallel to the yz plane, are very large compared to that of the edges of the web paralleled to the xy plane in Figure 3. The dimensions ofa typical web span (see Table 5) are 5m in length, 0.3048m (l ft) in width and 5.08 x 104 m (0.02 in) in thickness. The ratio of the crosssectional area of the yz plane to the crosssectional area of the xz plane is 600: I. Because there is very little heat transfer area on the edges, the temperature gradient toward the edges is insignificant, and the heat conduction in the y direction can 35  also be neglected. Three terms are left in the energy equation representing the forced convection due to the motion of the web in the machine direction, the heat conduction toward the large surfaces and the heat conduction in the machine direction due to a temperature gradient. Rewriting the heat fluxes in terms of temperature gradients yields c (3.30) A dimensionless group called the Pechlet number represents the ratio of the heat transfer by forced convection to that by conduction. The Pechlet number can be obtained from a dimensional analysis ofEq. (3.30) if the thermal conductivity is assumed constant. The steps of nondimensionalization of the energy equation can be found in Section 3.4. (3.31 ) The heat transfer is dominated by forced convection rather than conduction if the Pechlet number is far greater than I (Baird, 1995). Since most web systems have speeds on the order of a few hundred to a few thousand feet per minute, the Pechlet number typically ranges from a few hundred to tens of thousands. Therefore, the heat conduction in the machine direction is negligible. The final governing energy equation for the web system described in Section 3.1 and 3.2 emerges as: ~ or 8 ( or) pC v = k . p z 8z Ox Ox 36 (3.32) 3.3.4 Initial and Boundary Conditions The temperature at the beginning of the web span is specified as To. T=To @ z=O. (3.33) The temperature gradient and the heat flux in the x direction are symmetric about the centerline ofthe thickness (Assumption 11). x T~ __ ._. . . ..__. z b L Figure 5. A Web Span with Initial and Boundary Conditions As shown in Figme 5, the zaxis is defined along the centerline of the thickness of the web span to obtain a simpler boundary condition. The heat flux in the x direction is zero at the center line, which gives the first boundary condition. The second boundary condition is on the top or the bottom smface of the web where the web loses heat to the environment by the means of surface convection. This boundary condition is established by combining Fourier's law of heat conduction and Newton's law of cooling. where 37 (3.34) > h is the surface heat transfer coefficient, Tb is the temperature of the web at the bOlllldary exposed to air, Ta is the ambient air temperature. The initial and boundary conditions for the halfthickness control volume are summarized as, I.C. B.C. or =0 Bx @ z=O, @x=O, @x=b. (3.35) 3.4 Nondimensionalization The governing energy equation, Eq. (3.32), and its initial and boundary conditions, Eg. (3.35), need to be converted to dimensionless forms for general analysis, computation implementation, and future applications. A set of dimensionless variables and their derivatives are defined for these purpose. x z TT a= a b' ~ = L' e= T T ; o a aa=Bx az aT a~=L' 88= (3.36) b ' To Ta A few useful partial derivatives result from these definitions. 38 • These partial derivatives are substituted into Eq.(3.22) to obtain the dimensionless governing energy equation: (3.3 7) (3.38)  where blL is the ratio of web dimension, i.e. the half web thickness over the length of the web span. If the thermal conductivity is independent of temperature, k can be moved out [ PCpv zbJ ofthe partial derivative, and the PecWet number k :would appear on the lefthand side of the equation, together with the ratio of web dimension. Section 3.5 contains the derivation of an analytical solution of Eq. (3.38) when all material properties are asswned constant. Substituting the dimensionless variables and the derivatives to Eq. (3.35) gives the initial and boundary conditions in dimensionless form. I.e. B.C. 8=1 as =0 oa as hb =8 oa k @ p=O, @a=O, @a=l. (3.39) Another dimensionless group, the Nusselt number naturally appears. 39 (3.40) The dimensionless governing equation, shown in Eq. (3.38), is a second order nonlinear partial differential equation with nonlinear boundary conditions, Eq (3.40). The numerical solution methods and techniques are presented in the next chapter. 3.5 An Analytical Solution for Constant Properties If the material properties are assumed independent of temperature, Eq. (3.38) becomes (3.41 ) where the Pechlet number (Pe) and the web dimension ratio (Rc = b/L) are constants. The initial and boundary conditions remain in the same form as in Eq. (3.39), but the Nusselt number (Nu) is now constant. The assumption of constant properties has made the governing energy equation linear. An analytical solution is, therefore, possible for the energy equation under this assumption. Since the dimensionless temperature depends on two variables, a and P, the coupled dependency can be separated ifeis represented by the product of two new variables that are functions of only one ofthe independent variables. The method is known as separation of variables. For a parabolic partial differential equation, it has been proven that the solution can be expressed in this manner (Finlayson, 1980). Let the two new variables be x =X(a), and y = y(p). (3.42) Applying the separation of variables gives 40 e(a,p) =X(a)· y(p). The partial derivatives ofe are obtained by differentiation ofEq. (3.43). (3.43) and ae dX =y aa da ' 00 dY ap =XdP , (3.44) Substituting Eq. (3.44) into Eq. (3.41) to replace the partial derivatives ofe yields dY d2X (Pe .Rc)XdP = Y da2 . The dependent variables are separated. Another new function is defined to split Eq. (3.45) into two separate equations. I dY I d2X 2 (Pe·Rc)==A, . Y dP X da2 The solution of the first resulting ordinary different equation, dY A,2 = Y dP Pe·Rc' is very straight forward. The pdependent component ofe is where CI is the integration constant. The second resulting equation from the separation of variables is (3.45) (3.46) (3.47) (3.48) (3.49) Eq. (3.49) can be easily solved by examining its characteristic equation. The solution is 41 p x = C2 sin(Aa) +C3 cOs(Aa) , (3.50) c where C2 and C3 are integration constants. The first boundary condition is applied to solve for the integration constant. @a=O, Since Y =0 is a trivial solution, Be dX aa=O=Y. da (3.51) Carrying out the remaining arithmetic of Eq. (3.52) gives Thus, the first tenn in Eq. (3.50) vanishes and S can be expressed as: ).2 S =X(a)· y(p) =C4e PeRc~ cos(Aa), (3.52) (3.53) (3.54) where C4 = C1 C3• C4 and Aare the two remaining unknown variables in Eq. (3.54). The second boundary condition is introduced to solve for A. @a=l, Be au =Nu·S. (3.55) Substituting 8 in Eq. (3.55) with Eq. (3.54) gives Further simplification of Eq. (3.56) yields Nu tan(A)  ' A . Since tan('A) is a periodic function, there is more than one solution to Eq. (3.57), or 42 (3.56) (3.57) n = 0,1,2 .... (3.58) Due to the multiple solutions of A, Eq. (3.54) is transformed into a series: The last unknown variable C4n can now be solved by applying the initial condition: (3.59) Eq. (3.59) becomes @ p=O, 8=1. (3.60) (3.61) The series form of Eq. (3.61) prevents further reduction, unless the orthogonality of the cosine function is noticed. The orthogonality of the cosine function states: for m * n, (3.62) Therefore, Eq. (3.61) is multiplied by CO~Ama), then integrated. (3.63) According to Eq. (3.62), only when m = n, the right hand side ofEq. (3.63) is not zero. Further reduction can be done in the following steps: (3.64) 43 and, (3.65) • Substituting Eq. (3.65) back to Eq. (3.59) yields the analytical solution of the governing energy equation with constant material properties. where (3.66) (3.67) This analytical solution was used to verify the numerical solution in Chapter V. Parametric studies based on the dimensionless groups can also be conducted using Eq. (3.66) and (3.77). Eq. (3.41) has been solved analytically by Carslaw and Jaeger (1959). Carslaw and Jaeger used an infinite series approach for a more general initial condition: @ ~=O, (3.68) After proving the orthogonality of X functions, Carslaw and Jaeger took advantage of the following equation in their derivation. Their result has a different appearance: form;t:. n. (3.69) (3.70) When the thickness (L) in Eq. (3.70) is replaced with the dimensionless length, 1, and 44  sin(Am) f(a) with the dimensionless initial condition, 1, the integral becomes If the Am Nusselt number is substituted using Eg. (3.67), rearrangement ofEq. (3.70) results in Eq. (3.66). This confirms the validity of the analytical solution. 3.6 Summary In this chapter, the governing mass balance, Eq. (3.3), and force balance, Eq. (3.14), were developed for the single open span presented in Section 3.1. These two governing equations were combined to give Eq. (3.15), which can be used to calculate the web tension given the temperature profile ofthe web. The governing energy equation together with the boundary conditions was then developed, simplified and presented in the dimensionless form, Eq. (3.38) and Eq. (3.39). The solution ofEq. (3.38) gives the temperature profile of the web material under nonisothermal condition. However, these equations represent a highly nonlinear partial differential equation problem that can only be solved by numerical methods. Under the assumption of constant properties, an analytical solution represented by Eq. (3.66) and Eq. (3.67) can be found for the governing energy equation. Although such solution is not practical in the calculation of web tension because the web properties are indeed affected by temperature, the analytical solution provides a basis for the parametric study of the heat transfer phenomena in the single span web system. 45 • • CHAPTER IV NUMERICAL METHODS AND TECHIQUES An analytical solution to the dimensionless governing equation, Eq. (3.38), can not be accomplished without the assumption of constant material properties. Since such an assumption results in the loss of the physical significance of this study, the challenge was to select a reliable and convenient numerical method. Section 4.1 discusses the selection of the numerical method and the implementation procedures. This is followed by the formulation of Eq. (3.15) in Section 4.2 to accommodate the discrete temperature distribution in tension calculation. 4.1 Selection of Numerical Techniques The dimensionless governing equation, Eq. (3.38), is a highly nonlinear partial differential equation due to the temperature dependency of the material properties in the coefficients. The governing equation has two independent variables that originated from the two spatial variables in the transport and thickness directions. The second order derivative represents the heat conduction in the thickness direction. The first order term describes the forced convection in the transport direction. Since the boundary condition in the conduction (thickness) direction is in the form of the normal gradients on the 46 = boundary, and the boundary condition in the convection (transport) direction is a specified value at the beginning end, the governing partial differential equation can be treated as an initial value problem with one spatial independent variable. As discussed earlier in Chapter n, there are a number of numerical techniques for the solution ofthis nonlinear initial value, partial differential equation. The formulation of these methods for any nontrivial problem is a complicated process that is highly problem dependent. One minor adjustment in the model development phase can render the previous numerical formulation useless, and require new derivation and new techniques. Web Transport System (WTS) is a complicated simulator utilizing the idea of seven primitive elements and capable of dealing with the viscoelastic effects under either transient or steady state conditions. As research on other nonideal effects progresses, WTS will need to provide numerical solutions to individual or systems of ordinary and partial differential equations. Combining the integration computation into a few modules that are reliable and versatile, will greatly benefit the development process by saving tremendous time and effort in solving these separate problems. PDEONE developed by Madsen and Sincovec (l975a, I975b) emerges as the natural choice for initial value Partial Differential Equations (PDEs) because PDEONE uses the wellstudied method of lines (MOL) to provide centered differencing in the spatial domain. When interfaced with a robust integrator for ordinary differential equations (ODEs), PDEONE becomes a powerful integrator interface for initial value PDEs with a wide variety of initial and boundary conditions. 47 • b 4.1.1 Method of Lines and PDEONE The method of lines converts an initial value PDE or systems of initial value PDEs into a sets of coupled first order ODEs by finite difference approximations ofthe spatial derivatives. The implementation procedure is best illustrated with an example. A onedimensional rod is initially at a temperature, TB. The temperature at one end of the rod is changed to TL when t > 0, and the other end is kept a TB. The temperature variation ofthe entire rod over time is of interest. This example problem can be expressed mathematically by the following governing equation, initial and boundary conditions. aT a2T =D' (4.1) at ax2 ' I.e. T(x,O) = TB , (4.2a) B.C. T(O,t)=TB , (4.2b) T(L,t) = TL . (4.2c) where T is the target temperature function; x is the spatial variable; t is time; D is a diffusion coefficient; TB and TL are the initial and boundary values of the target function. The first step is to discretize the spatial domain, i.e. define the node points. For illustration purpose, five node points are used for this example. The initial and boundary conditions of this problem and the node points are shown in Figure 6. The distance between two adjacent node points is d. The second step is to discretize the partial differential equation at each interior node point. The second partial derivative with respect to x at each node is approximated 48 c  x ~~ L ~ Ti+1 d Ti I , i T I i1  "" '... t Figure 6. Node Point Layout of the MOL Example Problem in terms of the neighboring node values. The partial derivative with respect to time becomes a total derivative at each node point. Applying finite difference approximations at any interior point, i, yields the generalized ODE for the interior points. ( i = 2,3,4). (4.3) The boundary conditions for this problem are then used for the values ofthe exterior node points. If the boundary conditions are given in the form of normal derivative with respect to x, the ODEs for the nodes adjacent to the exterior nodes have to be modified to satisfy the boundary conditions. For the example problem, the MOL discretization yields the following system of ODEs. 49 (4.4a) (4.4b)  (4.4c) (4.4d) (4.4e) & The Method of Lines has generated three coupled ODEs, Eg. (4.4b), (4.4c) and (4.4d), which can be integrated simultaneously using the Gear method software, LSODE. The ODE integration method varies depending on the stability of the problem. A brief discussion on LSODE's solution ofODE systems can be found in Section 4.1.2. It is worth noting that the initial condition of the example problem, is now transformed to the discrete format. at t =o. (4.5) For more complex initial value PDEs than the example problem, a more general implementation scheme is necessary to account for the unlimited number of mathematical structures. According to Madsen and Sincovec (l975a), PDEONE is capable of such a task for the following structure that represents a wide class of physically realistic problems. 1 8 (c aT,) 1 a( c 8TNPDE )) 7 ax x Dk,l ax '''''7 Ox X Dk,NPDE Ox ' k = 1,2, ... , NPDE, (4.6) where NPDE denotes the number of initial value partial differential equations; and c is 0, 1 or 2 depending on whether the problem is in Cartesian, cylindrical, or spherical 50 j52 coordinates, respectively. The spatial domain is in the interval of (a, b), and the time domain is (to, 00). The initial conditions are & X E [a, b], k = 1, 2, ... , NPDE. (4.7) The boundary conditions are All ofthe coefficient functions (Uk, ~k, Yk, DkJ, fk, and ~k) must be at least piecewise continuous functions of an of their respective variables. Dkj (k, j = 1, 2, ... , may be functions oft and T, but only functions ofx othelWise. PDEONE can automatically handle all three types of boundary conditions, i.e. Dirichlet (~k = 0), Neumann (Uk = 0), or mixed (Uk:l: 0, ~k :I: 0). In addition, the initial condition need not be consistent with the boundary conditions. In other words, PDEONE does not require the initial condition to satisfy the boundary condition as x approaches either boundary a or b. The derivation ofthe finite difference approximation use by PDEONE is beyond the scope of this thesis, the details can be found in the paper written by Madsen and Sincovec (l975a). 4.1.2 LSODE Most integrators can typically solve initial value ODE systems of the form i = 1, 2, ... , NEQ. (4.9) 51 >  ODE integrators are often classified as being designed to solve either stiff or nonstiff equations. If there are more than one firstorder differential equation, the possibility of a stiff set of equations arises. A problem is considered stiff when the dependent variables change in two or more very different scales of the independent variable. Most nonstiffintegrators used the Adams predictorcorrector methods of up to twentieth order (Gear, 1971). The advantage of the integrators using the Adams method is that they usually do not require the solution of large matrix problems in their solution process since a fixed point iteration is used to solve the nonlinear equations. This remains true even when these nonstiff integrators are used with the PDEONE interface (Madsen and Sincovec, 1975a). The main disadvantage of the Adams method is that the stability regions for the higher order fonnulas are small and hence can often require very small time steps to maintain stability and accuracy. Because even simple PDE problems can become stiff as the number of spatial mesh points increases, nonstiff integrators are not always adequate. Stiff integrators usually use methods modified from Newton's method for solving the nonlinear equations (Madsen and Sincovec, 1975a). The iterative algorithm of these modified Newton's methods uses the Jacobian matrix Of / aT, and therefore repeated solutions of nontrivial matrix problems are required. The most widely used stiff integrators are based on Gear's stiffly stable fonnulas (Gear, 1971), which are generalizations of the backward difference methods. Due to the limitation of the nonstiff method and the computation requirement of the stiff method, an integrator with both nonstiff and stiff capabilities is necessary to provide flexibility for different applications. LSODE (Livennore Solver for Ordinary Differential Equations) is a powerful and 52   robust ODE integrator developed by Hindmarsh (1983) as an upgraded version of the older GEAR and GEARB packages. LSODE solves stiff and nonstiff systems represented by Eq. (4.9). In the stiff case, LSODE uses the Backward Differentiation Fonnula (BDF) methods treating the Jacobian matrix as either a fun or a banded matrix, and as either usersupplied or internally approximated by difference quotients. The resulting linear system is solved by direct methods (LV factor/solve). Adams predictor corrector methods are used for nonstiff cases without using the Jacobian. matrix. LSODE is a collection of subroutines that perfonn the following tasks: initialize integrator parameters, check for errors, set time integration method coefficients, perfonn a single step of the integration and maintain a user specified time integration accuracy, select the proper time step size and order in a dynamic way throughout the problem, and set up and solve any resulting linear algebraic systems. LSODE requires the user to provide the right hand side functions (f) ofEq. (4.9) and a main program to set up integration parameters and initiate the ODE integrator. For the initial value PDE problem derived in Chapter III, the PDEONE interface subroutine is used to convert the PDE to a system of ODE. The right hand side functions can be fonned and evaluated by PDEONE. The next section explains the details of interfacing LSODE, PDEONE, the initial value PDE, and the boundary conditions. 4.1.3 Overall Algorithm Structure and Subroutine Arguments The numerical techniques discussed in the previous section work in a straight task pattern as shown by a layout of the numerical approach in Figure 7. The initial value 53 I Initial Value PDE Eq. (3.38) & (3.39) ~ PDEONE Semidiscretization (MOL) I A System of Initial Value ODEs ~ LSODE ODE Integration (Adams/BDF) Numerical Solution Figure 7. Layout of the Numerical Approach POE and the boundary conditions, described by Eq. (3.38) and Eq. (3.39), are semidiscretized in the PDEONE subroutine using the method of Jines. The result ofthe POEONE semidiscretization is a system of initial value ODEs. The package ODE integrator, LSODE, solves this system of initial value ODEs by either the Adams predictor corrector methods or the Backward Difference Formula. The overall algorithm structure works in the reverse order of the task pattern due to initiation calls and output data processing. Figure 8 illustrates the structure of the overall algorithm. The main program initializes the problem by defining physical variables, setting parameters, calling the integration driver routine. The main program 54 also sends the numerical integration results to the tension calculation routine, and processes the output. Main Program Initialize, Call DRIVER, ...... Tension Calculation OutDut Results LSODE DRIVER ODE INTEGRATOR ROUTINES LSODE by Hindmarsh D Diffusion Coefficients MOL ROUTINES PDEONE, PDETWO by Madsen et at .E Partial Differential Equations Figure 8. Overall Algorithm Structure lBNDRY Boundary Conditions  The driver routine of LSODE makes calls to different subroutines of the LSODE package. The LSODE subroutines that require information of the system of ODEs will make subsequent calls to PDEONE. For each time step in the integration process, the PDEONE is called to perfonn semidiscretization in the spatial domain. In the actual interface, PDEONE is passed as a functional argument in the calls from the main program 55  to the LSODE driver, and from the LSODE driver to the LSODE subroutines. The LSODE subroutines recognize the first argument as the right hand side functions of the ODE system, fin Eq. (4.9). PDEONE requires the initial value PDE to be defined in three subroutines, D, F and BNDRY. Table 1 lists the arguments and their array size used in the communication between PDEONE and the three user defined subroutines. The first subroutine, D, evaluates the diffusion coefficients Dkj in Eq. (4.6). For heat transfer problems, the diffusion coefficient is the thennal conductivity. Argument Array Size Related Subroutines t Scalar D, F, BNDRY NPDE Scalar D,F,BNDRY x NPDE D,F,BNDRY U NPDE D,F,BNDRY DVAL NPDEbyNPDE D UX NPDE F DUXX NPDE F FVAL NPDE F ALPHA NPDE BNDRY BETA NPDE BNDRY GAMMA NPDE BNDRY Table 1. Arguments of the Three User Defined Subroutines 56 pas The temperature dependent function of web thermal conductivity needs to be enclosed in this subroutine. The input arguments of subroutine D consist of the number of POEs (NPOE), the independent variables (t and x), and the dependent variable (U). The subroutine returns the value of the diffusion coefficients (OVAL). The second subroutine, F, evaluates the right hand side functions fk in Eq. (4.6). As can be seen from Eq. (4.6), subroutine F requires first and second order derivatives (UX and DUXX) of the dependent variable in addition to the input arguments in Subroutine D. PDEONE approximates the first and second derivatives starting from the boundary conditions. Subroutine F only has to contain the algebraic fonnula in the format specified by Eq. (4.6). .The last user defined subroutine is BNDRY for boundary conditions evaluation. Subroutine BNDRY requires the same input arguments as subroutine D, and returns the values of the coefficient functions in Eq. (4.8), Uk, Pk, and ¥k. PDEONE makes calls to these three user defined subroutines to semidiscretize the POE problem and generate a system of ODEs. Therefore, whenever LSODE and its subroutines call PDEONE, the argument of PDEONE has to contain the number of spatial points (NPTS) and the vector variable (UOOT) of the size NPDE by NPTS to evaluate the right hand side of the ODE system. In addition, the main program call to the LSODE driver routine must include 17 arguments that are needed for problem definition, error control, method selection and optional inputs. Hindmarsh (1983) gave the detailed definition and application of these LSODE arguments. Table 2 contains brief definitions of these 17 variables and their respective values for the nonisothermal problems developed in Chapter III. 57 > Argument Definition Value Used f RHS ofODE system PDEONE NEQ Number ofODE NPDE by NPTS (1 x 101) y Dependent variable Dimensionless temperature (U) t Time Time tout Point ofdesired output Every 11100 of total length itol Tolerance method flag 1 for global absolute tolerance rtol Relative tolerance 1.0 x 105 atol Absolute tolerance 1.0 x 105 itask Task flag 1 for normal computation istate Calculation status flag 1 for normal initiation iopt Optional input flag ofor no optional input rwork Real working array Uninitiated when call. lrw Declare rwork length 11132 for NEQ = 101 iwork Integer work array Uninitiated when call liw Declare iwork length 121 for NEQ = 101 Jac User supplied Jacobian Not used mf Integration method flag 22 for BDF wi internal Jacobian Table 2. Arguments ofLSODE 58  4.2 Tension Calculation The tension level within the free span can be calculated from Eq. (3.15). Since the numerical solution of the governing energy equation generates a discrete temperature profile for the web material, the remaining task is the application of Eq. (3.15) to each small segment in the transport direction (z). Applying Eq. (3.15) to the ith segment yields (4.10) By definition, elongational strain in the transport direction for the ith segment is (4.11 ) where Lo is the length of the segment before elongation; El is the length of the segment after elongation. Web velocity and tension bear no subscript, because web velocity is constant under steady state conditions for an elastic material. Tension is also constant within a span due to Assumptions 4 and 9. The constant km can be eliminated by Eq. (3.3). Combining Eq. (4.10) and Eq. (4.11) with the elimination of km results in (4.12) L If n is the total number of segments in the transport direction, then nLo is the overaU length before elongation. Summing the length of each segment after elongation gives the overall length, EI. S9  > (4.13) Dividing nLo throughout Eq. (4.13), then subtracting 1 from both sides leads to the total elongational strain of the web span, Ez. (4.14) where Ez can be calculated from the roller speed difference using the following equation derived by Shin (1991). v2  VI E z =::...._' VI (4.15) ·Finally, with the simplification ofEq. (4.15), Eq. (4.14) can be rearranged for the tension within a span. Eq. (4.16) contains a summation of the thermal strain of each discrete element. According to the elastic assumption, these elements act as tiny springs connected in a series pattern. The overall tension of the connection is evidently contributed by the deformation or strain of each tiny spring. (4.16)  According to the discretization of the numerical approach shown in Figure 6, the nonisothennal problem has two dimensions, x and z. The z dimension (transport direction) is treated as time in the integration because Eq. (3.38) has the form of an initial value PDE. The result of the numerical integration is a temperature distribution in two dimensions. However, Eq. (4.16) needs only the temperature distribution in the z dimension due to the uniaxial tension assumption. In order to remedy this problem, an 60   average temperature over the distribution in the x dimension (thickness direction) is used in Eq. (4.16) to meet the uniaxial tension assumption. This approximation does not consider the shear stress created by different degrees of thermal expansion from the centerline to the surface of the web. 61 ..  CHAPTER V RESULTS AND DISCUSSION The nonisothermal model developed in Chapter III is capable of predicting the web tension based on the analysis ofthe temperature profile within a free web span. Since polypropylene (PP) is one of the most common web materials, PP was used to demonstrate the effect of nonisothermal conditions in web handling. Section 5.1 compares the dimensionless temperature results of the analytical solution of two test cases to the numerical results of the nonisothermal model. The comparison verifies the accuracy of the temperature profile computed by the numerical methods discussed in Chapter IV. Section 5.2 contains the temperature functions of the polypropylene properties obtained from published experimental data, as well as a list of other operating parameters. The tension module of the nonisothermal model is also verified under isothermal conditions. Section 5.3 presents the numerical simulation results of a physically realistic case under nonisothermal conditions. The resulting temperature profile of an open web span is examined, followed by a study of the thermal effects on web strains and tensions. The grid size of the numerical solution is optimized for computation time with 0.1% error tolerance in the nonisothermal tension prediction. Section 5.4 concludes this chapter with parametric studies of the model based on the dimensionless groups to provide insight to the nonisothermal effects on a web span. 62   5.1 Verification The web handling industry traditionally obtains the web temperature by online measurements, and controls tension using online empirical data. As with many aspects of web handling, there is a remarkable lack of data in the open literature to compare to the results of this study. In order to verify the numerical results of the nonisothermal model, two simple test cases were devised. The test cases use the following sets of constant parameters for polypropylene. The calculated values of the dimensionless groups for these parameters are also listed in Table 3 respectively. Case #1 represents a typical operating environment of a web system. Case #2 has the same parameters except that the web velocity is reduced almost 20 times to give an extremely low Pe in web handling. System Parameters of Test Case Case #1 Case #2 Density (p, kg/mJ ) 910 910 Heat Capacity (Cp , J/kgK) 2140 2140 Velocity (vz, ftlmin) 400 20.534 Thermal Conductivity (k, J/msK) 0.172 0.172 Surface Heat Transfer Coefficient (h, W/m2K) 100 100 Length of Web Span (L, m) 5 5 Half Thickness of Web (b. m) 2.54 x 10'4 2.54 x 104 Dimensionless Groups of Test Case Case #1 Case #2 Pechlet Number (Pe) 5844 300 Nusselt Number (Nu) 0.1477 0.1477 Size Ratio (Rc) 5.08 x 105 5.08 X 10'5 Table 3. Test Case System Parameters and Dimensionless Groups 63   The most important part of the nonisothennal model is the solution ofthe governing energy equation to obtain the temperature profile, which is then used to calculate the web tension. lfthe material properties are assumed constant, Eq. (3.38) becomes Eq. (3.41). Eq. (3.66) and Eq. (3.67) express one fonn ofthe analytical solutions. where . () 2 00 sm Am ~~ e= L: e Pe·Rc CO~A a) m=O Am + sin(2Am) m , 2 4 (3.66) (3.67) This analytical solution is in the form of an infinite series that converges quickly with two to three tenns for ~ values between 0 to 0.2. For larger Pvalues, the series is accurate up to the fourth significant figure using just the fust tenn. ~ S (0) e(I) e(2) e(3) e(4) LoS 0.000 1.0234 0.0286 0.0074 0.0033 0.0019 1.0000 0.025 1.0113 0.0122 0.0003 0.0000 0.0000 0.9994 0.050 0.9994 0.0052 0.0000 0.0000 0.0000 0.9943 0.075 0.9876 0.0022 0.0000 0.0000 0.0000 0.9854 0.100 0.9760 0.0009 0.0000 0.0000 0.0000 0.9751 0.125 0.9645 0.0004 0.0000 0.0000 0.0000 0.9641 0.150 0.9531 0.0002 0.0000 0.0000 0.0000 0.9530 0.175 0.9419 0.0001 0.0000 0.0000 0.0000 0.9419 0.200 0.9308 0.0000 0.0000 0.0000 0.0000 0.9308 0.225 0.9199 0.0000 0.0000 0.0000 0.0000 0.91'99 0.250 0.9090 0.0000 0.0000 0.0000 0.0000 0.9090 0.275 0.8983 0.0000 0.0000 0.0000 0.0000 0.8983 0.300 0.8877 0.0000 0.0000 0.0000 0.0000 0.8877 Table 4. Fast Convergence of the Analytical Solution 64  Table 4 shows the fast convergence of the analytical solution for the Test Case #1. The number in parentheses after 8 denotes the mlh term of the infinite series. The last column is the sum of all terms in the series. The results ofthe analytical solution for the test cases are plotted in Figure 9, along with the numerical results of the nonisothennal model. 1.0 Case # I e 09 .0.1.) E ~ 8. 0.8 e01) f Surface r.n r.n 01) 07 . "'2 0 'r;;  Analytical Solution c 01) E 06 • NIT Model 0 Pe = 5884 Nu =0.1477 05 00 0.1 02 0.3 0.4 0.5 06 07 08 09 1.0 Dimensionless Web Length (P) 1.0 0.9 Case #2 e 0.8 .0.1..) .2 0.7 .r.o.. 8. 0.6 . E  Analytical Solution 01) f 0.5 r.n • NIT Mode.1 r.n 01) 0.4 "'2 0 Pe = 300 Nu =0.1477 'r;; 0.3 . C 01) E 0.2 is 0.1 0.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Dimensionless Web Length (P) Figure 9. Dimensionless Temperature Profile of Test Cases 65  L The results of the analytical solution are shown by the solid lines to indicate the continuous nature of the analytical solution. The results of the numerical nonisothermaI model are represented by the square dots. As can be seen from Figure 9, the dimensionless temperature results obtained by both methods coincide with each other. For Case # I, the average relative difference between the numerical results and the analytical results is less then 0.01%. For the extreme Case #2, the average difference of the results is about 0.62% due to the almost zero dimensionless temperature for ~ greater than 0.6. In both cases, the dimensionless temperature results, on a scale from ato I, agree to the fourth significant figure, which is the specified error tolerance level of the numerical integration. This agreement indicates that the nonisothermal model can solve the governing energy equation accurately, and has been properly programmed for the test cases with constant properties. For a realistic heat transfer problem, the nonisothermal model incorporates the property functions in its numerical solution. The next section contains these property functions. 5.2 Property Functions and Parameters Four material properties that are functions oftemperature need to be supplied to the nonisothermal model before any analysis is possible. They are density, thermal conductivity, constant pressure heat capacity, and the Young's modulus. Predictive equations to estimate the temperature dependency of these properties can be found in Section 2.3. However, first hand experimental data and literature data are preferable. For 66  (5.1 )  L the purpose oftbis study, temperature functions of these four polypropylene properties are obtained by curve fitting the literature data. The unit of temperature is °C in all the following functions. The density data can be found in Mark (1996), and the function is P(kgl m3 ) = 8481.90x 1O2T3.05x 1O3T2 . The data for constant pressure heat capacity per unit volume was retrieved from Steere (1966), and the function is Cp(Call cm3 . K) =0.3397 + 1.028 x 1O3T + 1.164 x 1O5T2 . (5.2) Also retrieved from Steere (1966), the thennal conductivity function is k((all s· cm . K) =4.197 x 104 + 1.082 X 106T  1.045 x 108T2 + 1.090 x 1010T3 . (5.3) The Young's modulus function was regressed from stiffness data published by Reding (1958), and is E(psi) =93460  2385T + 34.86T2  0.2506T3 + 6.496 x 104 T4 . (5.4) The temperature function for the linear thennal expansion coefficient, pl, can be obtained from the density function using Eq. (3.12). Other operating parameters needed for the nonisothennal analysis include the dimensions of the web span, the roller velocities, overall surface heat transfer coefficient, initial and ambient temperatures. The operating parameters of a typical web span system are listed in Table 5. The surface heat transfer coefficient (h) was assumed a typical constant value of 100 W1m2K (Baird, 1995). The surface heat transfer coefficient is not required to be a constant by the numerical methods. As long as h is a continuous function of temperature 67 .......  Operating Parameters Value Unit Length of Web Span (L) 5 m Half Thickness of Web (b) 2.54 x 104 m Width of Web (w) 12 In Initial Web Crosssectional Area (Ao) 0.24 in2 First Roller Velocity (VI) 400 ftfmin Second Roller Velocity (V2) 401 ftfmin Initial Temperature (To) 120 °C Ambient Temperature (Ta) 25 °C Surface Heat Transfer Coefficient (h) 100 W/m2K Table 5. Operating Parameters of a Typical System VI = 400 [pm V2 = 401 fpm z ~ 11..4.t.  L =5 m.~I Side View x W=12in z .................._.__ ._.__ .  ..~ z y Top View Enlarged Side View L Figure 10. Web System Dimension Summary 68 or location, the user only needs to supply an additional subroutine to evaluate h in the same manner that the thermal conductivity (k) is evaluated for the boundary condition. Figure 10 summarizes the dimensions of the example web system. The system is a single span web with two driven rollers whose velocities are 400 ftlmin and 401 ftlmin, respectively. The initial web temperature (To) immediately after the contact zone ofthe first roller is at 120°C in all simulation unless noted otherwise. Similarly, the ambient temperature (Ta) is at 25°C unless noted otherwise. The polypropylene property functions in the temperature range from 25°C to 120°C are described by Eq. (5.1) through Eq. (5.4). The effect of the nonisothermal condition on web tension was the ultimate goal of this study. Since strain is directly related to tension, both strains and tensions under isothermal and nonisothermal conditions need to be examined. Under isothermal conditions, the web system does not undergo a temperature change, there is no thermal strain, and all physical properties remain constant. The only strain is the mechanical strain due to the different roller speeds. According to Eq. (4.15), the mechanical strain is 2.5 x 10'3 for the system defined by Figure 10 and Table 5. Consequently, Eq. (4.16) can be further simplified to (5.5) L When the ambient temperature (Ta) and the initial web temperature (To) are set to equal to each other, the nonisothermal model can be used to generate isothermal tension results that can be compared to Eq. (5.5). The results of the nonisothermal model integration and the results from Eq. (5.5) were plotted in Figure 11. The agreement of these results verifies that the tension module of the nonisotherrnal model has been 69 accurately programmed. Both Eg. (5.5) and Figure 11 show that the isothennal tension is only affected by the temperature via the Young's modulus. Figure 11 also demonstrates the sensitivity of web tension on temperature. 35 30 ~ 25 gc0 'Vi 20 cv r '"E 15 . v .c (5 ~ 10 5 0 20 30 • NIT Model  Equation (5.5) 40 50 60 70 80 90 100 Web Temperature/Ambient Temperature (0C) 110 120 Figure 11. Isothennal Tension at Different Temperature 5.3 Results Numerical simulations have been conducted for the system depicted by the property functions of polypropylene and the operating parameters presented in Section 5.2. Section 5.3.1 presents the temperature profile results. Section 5.3.2 contains the results and analyses on the tensions and strains at different initial and ambient temperatures. Section 5.3.3 examines the effect of grid size on the precision of results. 70 . 5.3.1 Temperature Profile Figure 12 is a plot of the dimensionless temperature profile for the nonisothermal case. The trend of the dimensionless temperature profile is similar to that of Test Case #1 shown in Figure 9. However, the dimensionless temperatures in Figure 12 decrease more gradually due to the incorporation of variable physical properties. At the surface end point of the web, the dimensionless temperature was 0.5928 in Test Case #1 with constant properties. The dimensionless temperature at the same location was 0.6508 in the nonisothermal model with variable properties. ''Q"") co "V; I:: Q) E i5 0.7 o a =0.0 " a = 0.2 )( a = 0.4 x a =0.6 • a = 0.8 + a = 1.0 0.\ 02 OJ 04 05 06 07 0.8 09 10 06 '~~~~~~~l 00 Dimensionless Web Length (P) Figure 12. Dimensionless Temperature Profile of a Single Polypropylene Web Span The difference may be explained by arguing that the material constant properties in the isothermal case tended to favor conduction over convection, hence, the web temperature profile was steeper in Test Case #1. Figure I in Chapter III shows that the 71  density of polypropylene increases slightly as temperature decreases while the heat capacity and thermal conductivity both decrease. Eq. (3.31) shows that the Pechlet number contains all three properties. (3.31 ) When the web cools, the decrease of k and the slow increase of p compete with the decrease of Cpo For polypropylene, the combined effect is a higher Pechlet nwnber indicating stronger convection at lower temperature in the nonisothermal case. Therefore, less heat is lost to the environment and the temperature at the end is higher in Figure 12. oX=o A X = b/5 )( X = 2b/5 ]I X = 3b/5 • X = 4b/5 110 ,, u t.... d..o.) a 100 .'".. do) 0. E do) f< 90 +X=b 80 '.~ o Web Length (m) Figure 13. Temperature Profile ofa Single Polypropylene Web Span The dimensionless temperature was converted back to actual temperature in 72 1.  r Figure 13. As can be seen in both Figures 12 and 13, the temperature profile over the length of web span appears almost completely linear, indicating a high convection to conduction ratio. This can be attributed to the high transport speed of the web leading to a relatively high Pechlet number. A high b/L ratio (Rc) will have similar effects, but the dimensions of the web span are usually not controllable operating parameters. Q • • • • • • • • • • • • • • • • • • • .. .. .. .. .. .. .. • • • .. .. .. .. .. .. .. • • x x x x x x x x x '" '" '" x x x '" '" J< J< + + + + + + + + + + + + + + + + + + + 0 c c c c c c c c c c c c c 0 c c c c 0 0 0 0 0 0 0 0 0 a 0 a 0 0 0 0 0 0 a • ~ = 0,0 .. ~ =0.2 J< ~ = 0.4 + ~ = 0.6 c ~ = 0.8 a ~ =1.0 l.l ~ 1.0 ~ .l.U.. ::I <..;.j. 0.9 lU C. E ~ 0.8 til til '!U C0 'c;; 0.7 l: !U E 6 0.6 0.5 0.0 0.2 0.4 0.6 0.8 1.0  Dimensionless Web Thickness (a) Figure 14. Dimensionless Temperature Profile of a PP Web Span over Thickness Figure 14 is a plot ofthe dimensionless temperature profile in the thickn.ess direction for the nonisothennal case. The temperature gradient over the thickness tends to increase gradually from the center to the surface of the web. This pattern is typical when Newton's cooling occurs at the boundary. Figure 14 also shows that the temperature gradient over the thickness of the web is fairly small because the actual measurement of web thickness is in the range of 0.01 in. Parametric studies of the 73    dimensionless temperature profile in the thickness direction, with the Pechlet and Nusselt numbers, can be fOWld in Appendix C. 5.3.2 Strains and Tensions Eq. (3.6) stays that the total strain (Ez) is the sum of the elastic strain (E) and the thennal strain (E1 ). The elastic strain is only caused by the tension in the web material. Since web material does not support compression, the elastic strain is always positive. On the other hand, because the linear thennal expansion coefficient is always positive, the thennal strain is negative if the web is cooling. Figure 15 is a plot of the strains over a range of initial temperature for the system defined by Figure lO and Table 5. 0008 T. = 25°C 0.006 • • Elastic Strain (E e ) • • 0.004 • • • •• • • •• • • • • • • • c:: 0.002 .(;j Total Strain (Ez) due to 1:: Mechancial Drawing r:rJ 0000 AA A • • • • • 0002 Thermal Strain (EI ) • • • ·0004 . 0.006 20 ]0 40 50 60 70 80 90 100 110 120 rnitial Web Tern peratu re, To (0C) Figure 15. Strains over a Range ofInitial Temperatures 74  As can be seen in Figure 15, the thermal strain becomes more negative if the initial temperature of the web is higher. As the difference between the initial temperature and the ambient temperature becomes greater, the elastic strain (ge), which is the quantity (Ez  f;t), becomes greater due to the decreasing (increasingly negative) value of gt. The web tension calculated from the elastic strain becomes larger also. The greater web tension is the result of the thermal contraction. From a mathematical perspective, since the total strain is fixed by the constant roller speeds, Eq. (3.6) requires a larger elastic strain to compensate for the decrease of thermal strain. In Figure 15, the total strain due to mechanical drawing is a constant value of 2.5 x 10'3 , which is the same as in the isothermal Test Case #1. Figure 15 also demonstrates the opposite directions of the elastic strain and thermal strain as the initial web temperature increases. The tension in the open span is also the sum of two components: tension due to mechanical deformation and the tension due to thermal shrinkage. For the web system defined by Figure 10 and Table 5, the nonisothermal model has performed the numerical integration over the initial temperatures (To) from 25°C to 120°C. The resulting tensions are plotted in Figure 16. When To = 25°C, the tension due to thermal strain is zero because this initial temperature is the same as the ambient temperature and there is no heat transfer. The tension due to mechanical strain and the overall nonisothermal tension have the same value of 31.15 lbf. As To increases, the tension due to thermal strain first increases, and reaches its maximum at about To = 110°C, then starts to decrease. The first increasing behavior results from the increasing thermal strain due to higher initial web temperature. 75  ...... T. = 25°C I •• • • • • • • .. • • • .. .. .. .. .. .. • • • • • • • • • Nonisothermal Tension • • Tension due to Thermal Strain • • • • • Elastic Tension  35 30 25 ~ :9 20 '" t:: 0 'Vi vt:: 15 flO o 20 30 40 50 60 70 80 90 100 110 120  Figure 16. Tensions over a Range ofInitial Temperatures Figure 16 also shows that the thennal strain becomes more negative. The negative sign of the thennal strain indicates the web is shrinking, At about To = 11 DOC, the Young's modulus has decreased to a point to offset the effect of thennal shrinkage. The tension due to thermal strain, which is a product of the Young's modulus and the thennal strain, peaks out and turns downward. The rapid decreasing of the Young's modulus with temperature is also reflected by the fast decreasing trend ofthe tension due to mechanical strain. The nonisothennal tension, which is the sum of the tensions due to mechanical strain and thennal strain, first decreases with To, then starts to increase slowly at about 50°C, finally turns to decrease at about 9DoC. This pattern indicates the competing effects of the thennal strain and the Young's modulus under the nonisothennal condition. In this case study, the thermal strain increases and the Young's modulus decreases with 76  higher initial temperature. The two competing effects tend to lessen each other's impact on the nonisothermal tension at the mid temperature range ofTo, i.e. 40°C  100°C. The nonisothermal tension is expected to decrease drastically if the initial temperature is beyond 110°C, because the Young's modulus is so small that both the tensions due to mechanical strain and due to thermal strain are declining. 20 25 30 ,15 40 To = 120°C • • • • • · • • • • • • • Nonisothermal Tension .  • Tension due to Thermal Strain • Elastic Tension o 15 10 30 § 15 .Vj ~ 20 25 Figure 17. Tensions over a Range of Ambient Temperatures The effect of the ambient temperature, Ta, on the tensions has also been studied. Under the condition defined by Figure 10 and Table 5, the simulation tension results were plotted against ambient temperatures in Figure 17. The nonisothermal tension and its component tensions are less sensitive to Ta than to the initial web temperature (To), because the change of Ta is rather small when compared to the already large difference between To (= 120°C) and Ta. In other words, the change in the driving force of heat  77   transfer is small. The decreasing trends ofthe nonisothermal tension and both the components with increasing Ta can be interpreted in terms ofthe two most significant effects of the nonisothermal condition. When the level of the web temperature profile rises slightly with higher Ta., the Young's modulus of the web (E) becomes slightly smaner. Therefore, Figure17 shows that the tension due to mechanical strain decreases slowly with Ta. On the other hand, because the initial temperature difference, (To  Ta), becomes smaller with higher Ta, the temperature profile becomes flatter. Since the thermal strain indicates how far the web temperature decreases from the initial temperature, the flatter the temperature profile, the smaller the thermal strain. Consequently, as a product ofE and Et,the tension due to thermal strain decreases with increasing Ta. Instead of competing, E and Et reinforce each other when the ambient temperature rises. In Figure 17, the nonisothermaJ tension shows a faster decrease because it is the sum of the decreasing tensions due to both mechanical strain and thermal strain. 5.3.3 Grid Size In this section, the grid size used in the numerical integration and the tension calculation is examined to achieve an acceptable error tolerance. Since different web systems will have different physical properties and different error requirement, the web system and the error tolerance must be defined before adjusting the grid size. For the purpose of this study, the web system for the grid size analysis is again described by Figure 10 and Table 5. 78 The grid size (or mesh size) in this discussion is in the fonn ofNPTS x NSTEP. NPTS is the number of points over one half ofthe web thickness because the numerical integration is only perfonned from the center to the surface of the web in the thickness direction. NSTEP is the number of grip points in the transport direction. Since the transport direction is considered as time by the integration routine, LSODE, NSTEP is originated from the number of time steps. However, the actual integration steps are internally generated by PDEONE to match the specified error tolerance, which was set at an absolute value of 1.0 x 10.5 to give 4 significant figures accuracy for the dimensionless temperature. Therefore, NSTEP is only used to control the points where the integration routine will output the dimensionless temperature. The dimensionless temperature at the surface end point converged to 0.6508 for 20 meshes, i.e. NPTS = 21. NSTEP was only set at 2 for a single output point. The temperature accuracy was controlled internally by the integration algorithm. This demonstrates that the PDEONELSODE numerical integration approach is a stable and converges rapidly. However, since the tension calculation involves the summation ofthe properties at each grid point, the mesh size needs to be large enough to account for the sensitive temperature function of the Young's modulus. The nonisothermal tension of the web system defined by Figure 10 and Table 5 was solved repeatedly for different mesh sizes. For To = l20De and Ta = 25°C, the nonisothennal tension converged to 24.64 lbf. The mesh size for this value was found at 501 x 5001. Further increase of the mesh size shows no change in the nonisothennal tension up to the fourth significant figure. Therefore, 24.64 lbf is used as the reference tension to obtain the percent error of the 79  other mesh sizes. The percent error is rather a precision error that an accuracy error. F 24.64 %Error = x100% 24.64 (5.6) 06 c .2 U ~ 0.4 ... 0.. c0 'r;; 0.2 C Cl.l f< lil E 0.0 Cl.l 5 0en ·c 0.2 0Z .5... 0 t: 0.4 W '#. • • • • • • 0 • 0 " " " " " " 0 • • • •. 0 0 Q • • NPTS = II • o NPTS = 21 .. .. NPTS = 51 o NPTS = 101 + NPTS = 201 o NPTS = SOl 0.6 100 1000 10000 100000  Number of Transport Direction Grid Points (NSTEP) Figure 18. The Nonisotherrnal Tension Prediction Error of Different Grid Sizes Figure 18 is a plot of the percent error in the nonisotherrnal tension for different mesh sizes. From Figure 18, less than 0.1% error in the nonisotherrnal tension prediction can be achieved with a grid size of 51 x 501 or larger. If 0.1% error in precision is required, 51 x 501 is the recommended mesh size. The series ofNPTS = 101 (represented by the circles), and the series ofNPTS = 201 (represented by the crosses) coincide almost everywhere, except at NSTEP = SOland 2001. This slight fluctuation of the errors indicates how sensitive the tension is to the Young's modulus and the thermal strain. 80 ......  Since the nonisothermal tension is directly affected by the thermal strain, the percent error in thermal strain prediction can be used as another precision index. The reference strain is 4.053 x 103 , which is again obtained with a mesh size of 50 I x 5001. The definition of the percent error in thermal strain is similar to Eq. (5.6). • NPTS = II c NPTS = 21 '" NPTS = 51 o NPTS = 101 . + NPTS = 201 ., NPTS = 501 c • ~ • D • • • • '" D D ~ '" D D 0 0 ~ •'" 0'" 0'" 0'" 14 C I 2 .9 '0 :a~... 1.0 0.. c '(..;.i (j) 0.8 «i E ~ 0.6 .c: f ..!.::. 04 0 t: w ~ 0.2 0 0.0 סס 100 1000 10000 10 oo  Number of Transport Direction Grid Points (NSTEP) Figure 19. The Thermal Expansion Prediction Error of Different Grid Sizes Figure 19 shows that more grid points are necessary to keep the error less than 0.1 %. Consequently, the number of the grid points in the form ofNPTS x NSTEP is recommended to be at least 101 x 1001. The increased mesh size requirement for the same precision in thermal strain can be attributed to the competing effects of thermal strain and Young's modulus. It was concluded earlier from Figure 14 that these competing effects decrease each other's impact on the nonisothermal tension. 81 ......  5.4 Parametric Studies The nonisothermal model contains more than a dozen operating parameters and property functions. Individual parametric study only offers disconnected information on the general behavior ofthe nonisothermal model. Dimensionless groups of these operating parameters can better illustrate the effects of the parameters under nonisothermal conditions. During the derivation of the governing energy equation, two dimensionless groups, Pe and Nu, were encountered.  and hb Nu=k . (3.31) (3.40) b The Pechlet number (Pe) represents the ratio of the heat transfer by forced convection to conduction. In the nonisothermal model, the forced convection refers to the energy transported into and out of the span via the movement of the web. A portion of the heat brought in the web span control volume by forced convection is transferred to the web surface by conduction and lost to the environment by surface convection. The remaining portion of the heat is transported out of the control volume by the forced convection. According to Eq. (3.31), Pe is the gross representation of density, heat capacity, thermal conductivity, web velocity and thickness. The second dimensionless group is the Nusselt number (Nu) which shows the dimensionless temperature gradient averaged over the heat transfer surface. This definition is the same as the common definition ofNu, which is the ratio of surface convection to conduction. Both definitions can be visualized by combining Eq. (3.10) 82 ....  and the boundary condition on the surface as depicted by Eq. (3.39) (5.7) Similarly, the Nusselt number is the gross representation of surface heat transfer coefficient, web thickness and thennal conductivity. Since Pe and the Nu are dimensionless quantities, the dimensionless temperature is a natural choice for the parametric study. Figure 20 shows how the dimensionless temperature at the web center (ex = 0) corresponds to different Pechlet numbers. The Nusselt number used to obtain the data points on this plot is 0.1477, which was also used in Test Case #1. 09 10 Pe=8000 0.7 08 a = 0, Nu =0.1477 04 0.5 0.6 + + + 0.3 Pe=IOOO 0.1 0.2 Pe=300 Pe=50000 +•J•o.•. • I I ·li·S· · · · · · · · · · · · · · · · · · · · · · · · · •••••• .. •• Itt·.. ~ ... . + .. .. .. • • • : : • • • • • • l)e= 10000 + .. • • & .. • • • • I & ... • • • • A .. ... ..... ••••••• + • .. • • • • .. .. .. ... ... & • • + .. .. . Pe=5844· • • • .. ...... ... ......~ + ... & ... ;' • • • • + Pe=3000" .......... + ....... + ..... ... .to .... + + ... ... .. ... 1.0 CD '' eu 08 l aos l eu Co E 06 eu fo Cfl ~ t: 04 'r;; Ceu E 6 02 00 0.0 Dimensionless Web Length (13) Figure 20. The Effect ofPechlet Number on the Dimensionless Temperature A very large Pechlet number would correspond the web moving through the open  83 ......  span so fast that there is little or no time for conduction to the surface. The web temperature change is very small over the length for large Pechlet numbers due to the slow conduction. Figure 20 confIrms that the dimensionless temperature at the web center decreases slowly, and the profile is almost completely linear for Pe > 3000. For small Pechlet numbers, conduction is rapid compared to transport through the span. As shown in Figure 20, the smaller the Pechlet number, the faster the initial drop in dimensionless temperature. ln addition, the rate at which the dimensionless temperature decreases, slows at large dimensionless web length (P) because the conduction heat transfer driving force (T(b=O)  T(b=l) ) is reduced. Figure 21 shows how the dimensionless temperature profile at the web center (a = 0) is affected by the Nusselt number. The Pechlet number used to obtain the data points on this plot is 5844, which was also used in Test Case #1. 08 09 10 a = 0, Pe =5844 os 0.6 07 c c cacooooQCCnn~nnnn • • • • •• • • • • • • • • • • • • • • • • Nu=1 c c o • 0 0 0 • 0 • 0 0 • 0 c 0 • 00 0 • 0 • 0 0 • c • c Nu=O.OI 1 0 ~ I I I I I • I : : : • • • • • • • • • • • • • • • • • • • • • • • • • • • • o : 0 0 t t ~ •x •• •: : •••• •A • • • • • • • • • • • • • • • Nu=O.05 • 0 • • • z • • A A • A • • • • • • • • • 0. 0 0 •••• "'''' •••••••• Nu=O.1 • 0 0 • Z K * • * • • A • A • • • •• 0 00 Nu=O.1477 ••••••• • J( :I: K 02 0 0 Nu=5 c c c c c 0.0 0.0 01 0.2 OJ 04 e~ 0.8 ~ ll) 0 E 06 ~ en en ~ c: 04 'v; c: ll) E is Dimensionless Web Length (~) Figure 21. The Effect of Nusselt Number on the Dimensionless Temperature 84  ......  l Since the Nusselt nwnber is the ratio of surface convection to conduction, a very small Nusselt number represents that the convection on the surface is so slow that the web is almost insulated. The web temperature change is very small due to the slow convection. The result of a very small Nusselt number is similar to a very large Pechlet number in that little heat transfer occurs between the environment and the web. Figure 21 shows that the Nusselt number affects the dimensionless temperature in much the same pattern as the Pechlet number, except in an opposite manner. For a large Nusselt number, i.e. Nu> 0.5, the strong surface convection over conduction results in nonlinear temperature profiles. The quick initial drop in temperature resembles the effect of a Pechlet number less than 3000. The temperature profile analysis with the Pechlet number and the Nusselt number reveals that the heat transfer process in the web system involves three steps, forced convection, conduction and surface convection. Other than the temperature differences, each step is affected by some operating parameters and physical properties. Since the physical properties are in turn affected by the temperature, it is necessary to have a nonisothermal model capable to deal with this coupling relationship. 85 ......   CHAPTER VI CONCLUSIONS AND RECOMMENDATIONS 6.1 Conclusions This study was successful in modeling the effects of nonisothermal conditions on a single web span. The development of the nonisothennal model and numerical simulation provided better understanding of the heat transfer between the web and the environment. The significance of this work rests on the ability to predict the web temperature profile, hence, the local physical properties critical to tension control. Web handling systems that are under nonisothennal condition can use the nonisothermal model to accurately analyze the thermal effects, which the currently existing isothermal models and previous attempts can not capture. The major conclusions are summarized below. 1. A nonisothennal model was developed to account for the thermal effects in web handling systems whose temperature is different from the temperature of the environment. The temperature dependence of the needed material properties were incorporated into the nonisothennal model, which was verified in a test case using constant material properties. The resulting partial differential equation was solved using widely accepted and robust numerical routines, PDEONE and LSODE. 86 r  2. A modeling framework was established to accurately evaluate the temperature distribution within the web. Since the energy equation was solved independently, the model can be adjusted to solve problems with other nonideal effects in conjunction with other models, such as Guan's Viscoelastic Model (VEM). 3. Web tension was found to be significantly affected under nonisothennal conditions in two ways, by inducing a thermal strain and altering the physical properties. The larger the temperature difference between the web and the environment, the greater the predicted isothermal tension deviation from the nonisothermal tension. 4. Among the physical properties encountered in this study, the Young's modulus was the most sensitive to temperature. Since tension is directly proportional to the Young's modulus in an elastic model, accurate determination of the Young's modulus throughout the web is imperative. 5. The mesh size in the numerical simulation is as important as an accurate temperature function of the Young's modulus. Significant error in tension and strain prediction can arise if the mesh size is not large enough. Trial runs with different mesh sizes are recommended to meet specific precision requirement. 6. In a cooling process, the thermal strain is negative and the web is shrinking. The larger overall strain makes the nonisothermal tension larger than the mechanical tension calculated from the velocity difference of the rollers. In a heating process, the thermal strain is positive and the web is expanding. In order to avoid sagging, the mechanical strain has to be greater than the thermal strain to retain positive tension. 7. Two dimensionless groups, Pe and Nu, can be used to study the effect of operating parameters in the nonisothennal model over the temperature distribution. The 87 .... thermal effects on the web not only depend on the material properties and the temperature difference, but can also be controlled by the operating parameters. 6.2 Recommendations In the use of the nonisothermal model to predict web tension, it is recommended that the temperature functions of the properties be determined by experimental measurement, especially, that ofthe Young's modulus. Caution should also be exercised in selecting the number of grid points to ensure precise prediction. Future investigation should be aimed at reducing the limitations of the model posed by the following assumptions. Assumption 2: The web material behaves purely elastically. Assumption 3: The web system is operating at steady state. Assumption 10: Viscous dissipation is insignificant compared to surfaced convection. The VEM model developed by Guan (1995) is capable of solving isothermal, mUltispan, viscoelastic problems under both steady state and transient conditions. The merge of the VEM model and the nonisothermal model will eliminate Assumptions 2 and 3. The resulting model can have the combined benefit to solve web handling problems with both thermal and viscoelastic effects. However, there are several changes that have to be made in the model development phase. 1, The significance of viscous dissipation (Assumption 10) needs to be studied. 2, A time dependent term should be included in the energy balance for unsteady 88  state, leading to a partial differential equation with two spatial variables and one time variable. 3. The force balance in the VEM model should be redeveloped to account for thennal strain. 4. Simplifications made under the assumption of isothermal conditions in the derivation of the VEM model should be reworked. 5. The temperature dependence of the material properties should be reflected in the VEM model. 89   REFERENCES 1. Aklonis, J. J., MacKnight, W. J., and Shen, M., 1972, Introduction to Polymer Viscoelasticity, WileyInterscience, New York. 2. Ames, W. Fo, 1965, Nonlinear Partial Differential Equations in Engineering, Academic Press, New York. 3. Back, E. L., and Andersson, L. 1., 1993, "The Effect of Temperature on Wet Web Strength Properties," Tappi, Vol. 76, No.5, pp. 164172. 4. Baird, D. Go, and Collias, D. 1., )995, Polymer Processing Principles and Design, ButterworthHeinemann Series in Chemical Engineering, Boston. 5. Balasubramaniurn, S., 1997, "Development of Interactive Software to Study the Effects of Viscoelasticity in Web Lines," Master's Thesis, Oklahoma State University. 6. Bicerano, 1., 1993, Prediction of Polymer Properties, Marcel Dekker, Inc., New York. 7. Bird, Ro 8., Stewart, W. E., and Lightfoot, E. N., 1960, Transport Phenomena, John Wiley & Sons, New York. 8. Bird, R. B., Annstrong, R. c., and Hassager, 0.,1987, Dynamics of Polymeric Liquids, Volume 1, Fluid Mechanics, Second Edition, John Wiley & Sons, New York. 9. Boresi, A. P., 1965, Elasticity in Engineering Mechanics, PrenticeHall, Inc., Englewood Cliffs, New Jersey 10. Boyd, 1. Po, 1989, Chebyshev and Fourier Spectral Methods, SpringerVerlag, New York. 11. Brandenburg, G., 1976, "New Mathematical Models for Web Tension and Register Error," Proc. 1, 3rd International IFAC Conference on the Instrumentation and Automation in the Paper, Rubber and Plastics Industries, Brussels, May 2426. 12. Brezinski, J. P., 1956, "The Creep Properties of Paper," Tappi, Vol. 39, No.2, pp. 116128.  90 1. _   13. Burnett, D. S., 1987, Finite Element Analysis: From Concepts to Applications, AddisonWesley, Reading, Massachusetts. 14. Byrd, Y. L., 1972, "Effect of Relative Humidity Changes During Creep on Handsheet Paper Properties," Tappi, Vol. 55, No.2, pp. 247252. 15. Campbell, D. P., 1958, Process Dynamics, John Wiley & Sons, New York. 16. Canuto, C., Hussaini, M. Y., Quarteroni, A., and Zang, T. A, 1988, Spectral Methods in Fluid Dynamics, SpringerVerlag, New York. 17. Carslaw, H. S., and Jaeger, J. C., 1959, Conduction of Heat in Solids, Second Edition, Oxford University Press, Oxford. 18. Chandrupatla, T. R., and Belegundu, A D., 1991, Introduction to Finite Elements in Engineering, Prentice Hall, Upper Saddle River, New Jersey. 19. Dugdale, D. S., 1968, Elements of Elasticity, Pergamon Press Ltd., Oxford 20. Dunn, H. S., 1969, "Equivalent Circuit Representation of Web Conveyance Systems," Journal of Applied Mechanics, June, pp. 316318. 21. Ferry, J. D., 1970, Viscoelastic Properties of Polymers, Second Edition, John Wiley & Sons, Inc., New York. 22. Finlayson, B. A., 1980, Nonlinear Analysis in Chemical Engineering, McGrawHill, New York. 23. Forsythe G. E., and Wasow, W.R., FiniteDifference Methods for Partial Differential Equations, John Wiley & Sons, New York. 24. Gao, 1., and Weiner, 1. H., 1992, "Computer Simulation of Viscoelasticity in Polymer Melts," Macromolecules, No. 25, pp. 13481356. 25. Gear, C. W., 1971, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, Englewood Cliff, New Jersey. 26. Gehlbach, L. S., Kedl, D. M., and Good, J. K., 1989, "Predicting Shear Wrinkles in Web Spans," Tappi Journal, August, pp. 129134. 27. Gottlieb, D. and Orszag, S. A., 1977, Numerical Analysis of Spectral Methods: Theory and applications, S.lA.M., Philadelphia. 28. Grenfell, K. P., 1963, "Tension Control on PaperMaking and Converting Machinery," Proc., IEEE Ninth Annual Conference on Electrical Engineering in 91  p:s the Pulp and Paper Industry, Boston Massachusetts, June 2021. 29. Guan, x., 1995, "Modeling of Viscoelastic Effects on Web Handling System Behavior," Doctoral Thesis, Oklahoma State University. 30. Harper, C. A, 1975, Handbook of Plastics and Elastomers, McGrawHill, Inc., New York. 31. Hartnet, J. P., 1992, "Viscoelastic fluids: A New Challenge in Heat Transfer," Transactions of the ASME, Vol. 114, May, pp. 296303. 32. Hauptmann, E. G., and Cutshall, K. A, 1977, "Dynamic Mechanical Properties of Wet Paper Webs," Tappi, Vol. 60, No. 10, pp. 106108. 33. Henderson, G. M., Hung, 1. C. and Googe, 1. M., 1979, "Automatic Tension Control for Winding Superconducting Coils," Proc. Southeast Conf. Reg. 3 Conf., Roanoke, Virginia. 34. Hilado, C. 1., 1969, Flammability Handbook for Plastics, Technomic Publications, Stamford, Connecticut. 35. Hindmarsh, A C., 1972, GEAR: Ordinary Differential Equation System Solver, Rep. UCID30001, Rev.2, Lawrence Livennore Lab., Livennore, California. 36. Hindmarsh, A C., 1973, GEARB: Solution of Ordinary Differential Equations Having Banded Jacobian, Rep. UCID30059, Lawrence Livermore Lab., Livennore, California. 37. Hindmarsh, A C. 1983, "ODEPACK, A Systematized Collections of ODE Solvers," In Scientific Computing, Stepleman, R. S., et al. (eds.), NorthHolland, Amsterdam, pp. 5564. 38. Hindmarsh, A c., 1987, "Brief Description ofODEPACK  A Systematized Collection ofODE Solvers," http://www.netlib.org/odepack. 39. Hopfinger, A J., Koehler, M. G., and Pearlstein, R. A, 1988, "Molecular Modeling of Polymers. IV. Estimation of Glass Transition Temperatures," Joumal of Polymer Science: Part B 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


