

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


DEVELOPMENT OF NONLINEAR FINITE ELEMENT CODES FOR STABILITY STUDIES OF WEB MATERIAL WITH STRAIN STATE DEPENDENT PROPERTIES By HASAN H.YURTCU Bachelor of Mathematical Engineering Yildiz Technical University Istanbul, Turkiye 2003 Master of Science Mathematical Engineering Yildiz Technical University Istanbul, Turkiye 2006 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of DOCTOR OF PHILOSOPHY December, 2011 ii DEVELOPMENT OF NONLINEAR FINITE ELEMENT CODES FOR STABILITY STUDIES OF WEB MATERIAL WITH STRAIN STATE DEPENDENT PROPERTIES Dissertation Approved: Dr. J.K. Good Dissertation Adviser Dr. R.P.Singh Dr .A.K.Kalkan Dr.B.Binegar Outside Committee Member Dr.Sheryl A.Tucker Dean of the Graduate College iii TABLE OF CONTENTS Page LIST OF SYMBOLS .................................................................................................................. v LIST OF FIGURES .................................................................................................................. vii CHAPTER Page 1. CHAPTER I ............................................................................................................. 1 1.1 INTRODUCTION ................................................................................................... 1 2. CHAPTER II ........................................................................................................... 7 2.1 LITERATURE REVIEW ........................................................................................ 7 2.2 Theory of Membrane Instability .............................................................................. 7 2.3 Applications of Instability Analyses to Webs ....................................................... 11 2.4 Summary ................................................................................................................ 22 2.5 Statement of Proposed Research ........................................................................... 23 3. CHAPTER III ........................................................................................................ 24 3.1 FINITE ELEMENT EQUATIONS ....................................................................... 24 3.2 Two Dimensional Four Node Quadrilateral Elements .......................................... 24 3.3 Numerical Integration by Gaussian Approach ...................................................... 33 3.4 The Global Stiffness Matrix and Matrix Banding ................................................. 37 3.5 A [D] Matrix for Plane StressIsotropic Material and Partly Wrinkled Membrane Elements ................................................................................................................ 42 4. CHAPTER IV ........................................................................................................ 50 4.1 MISALIGNED ROLLERS .................................................................................... 50 4.2 Beisel’s Method for Modeling Wrinkles Due to the Misaligned Rollers .............. 50 4.3 A New Algorithm for Predicting Wrinkles Due to the Misaligned Rollers .......... 54 4.4 Comparison of Results with Previous Works ........................................................ 66 4.5 Improvement of Execution Time ........................................................................... 68 4.6 Automating the Code ............................................................................................. 72 4.7 A New Slack Edge Criteria for a Misaligned Roller ............................................. 84 5. CHAPTER V ......................................................................................................... 91 5.1 TAPERED ROLLERS .......................................................................................... 91 5.2 Description of Tapered Rollers and Beisel’s Method for Modeling Wrinkles Due to iv CHAPTER Page the Tapered Rollers ................................................................................................ 91 5.3 A New Algorithm for Predicting Wrinkles Due to the Tapered Rollers ............... 95 5.4 Automating the Code ........................................................................................... 101 5.5 Comparison of Results with Previous Works for the Tapered Roller ................. 103 6. CHAPTER VI ...................................................................................................... 109 6.1 HOLES IN WEBS ............................................................................................... 109 6.2 Mallya’s Method for Modeling Wrinkles Due to Circular Discontinuity ........... 109 6.3 A New Algorithm for Predicting Wrinkles Due to a Circular Discontinuity ...... 113 6.4 Automating the Code ........................................................................................... 119 6.5 Comparison of Results with Previous Works for Wrinkles Due to Holes .......... 121 7. CHAPTER VII..................................................................................................... 128 7.1 NON UNIFORMITIES IN WEBS ...................................................................... 128 7.2 Modeling Non uniformities with Commercial Finite Element Code COSMOS . 129 7.3 A New Algorithm for Predicting Wrinkles Due to the Circular Non uniformities132 7.4 Comparison of COSMOS and Excel VBA Results ............................................. 138 7.5 Automating the Excel VBA Code ....................................................................... 142 8. CHAPTER VIII ................................................................................................... 145 8.1 CONVERGENCE CHECK ................................................................................. 145 8.2 The Misaligned Roller Case ................................................................................ 146 8.3 The Tapered Roller Case ..................................................................................... 154 8.4 The Central Circular Void Case .......................................................................... 161 8.5 The Central Circular Non Uniform Web Case .................................................... 166 9. CHAPTER IV ...................................................................................................... 169 9.1 SUMMARY ......................................................................................................... 169 9.2 RECOMMENDATIONS FOR THE FUTURE STUDY .................................... 170 REFERENCES ....................................................................................................................... 172 APPENDIX EXCEL VBA CODE FOR MISALIGNED ROLLERS .................................... 178 v LIST OF SYMBOLS a Web Span Length E Young’s Modulus FE Finite Elements Fig. Figure MD Machine Direction CMD Cross Machine Direction R Radius of the roller t Web Thickness Tw Web line tension X,Y,Z Coordinates σ Stress ycr σ Critical Buckling Stress md σ Machine Direction Stress cmd σ Cross Machine Direction Stress ν Poison’s Ratio ε Strain μ Coefficient of Friction I, Iz Area Moment of Inertia mcr Critical Roller Taper T Web Tension Tslack Tension Where Slack Edge Occurs vi V, Vavg Velocity, Average Velocity ω Angular Velocity σ1, σ2 Principal Normal Stresses Ffriction Friction Force FE Finite Elements x,y,z Coordinate Directions σy Normal Stress εx, εy Normal Strains ε1, ε2 Principal strains vii LIST OF FIGURES Figure Page Figure 1. 1 Troughs and Wrinkles .............................................................................................. 2 Figure 1. 2 Regime I and Regime II Wrinkles ........................................................................... 3 Figure 2. 1 Boundary Conditions for Misaligned Roller .......................................................... 11 Figure 2. 2 Webb’s Model ........................................................................................................ 15 Figure 2. 3 Webb’s Model: Boundary Conditions and Applied Loads .................................... 16 Figure 3.1 Four Node Quadrilateral Element .......................................................................... 26 Figure 3. 2 Quadrilateral Elements Inξ ,η Space (Master Element) ....................................... 27 Figure 3. 3 Two Points Gauss Quadrature using 2x2 rule. ....................................................... 36 Figure 3. 4 A Simple Model for Global Stiffness Matrix ......................................................... 37 Figure 3. 5 A Quadrilateral Element and Nodal Displacements .............................................. 38 Figure 3. 6 Mohr’s Circle for Plane Strain ............................................................................... 44 Figure 3. 7 Behaviors of Wrinkling Membrane Elements [12] ................................................ 48 Figure 4.1 Beisel’s FE Wrinkle Model for Misaligned Roller ................................................. 51 Figure 4.2 L = 6’’ Span Results ................................................................................................ 53 Figure 4.3 L = 18’’ Span Results .............................................................................................. 53 Figure 4.4 L = 30’’ Span Results .............................................................................................. 54 Figure 4.5 The System That is Modeled................................................................................... 58 Figure 4.6 Misaligned Roller FE Wrinkle Model .................................................................... 59 Figure 4.7 Flow Chart for New Algorithm ............................................................................... 61 Figure 4.8 Convergences of P and Q ........................................................................................ 63 Figure 4.9 Excel Input for the Program .................................................................................... 64 Figure 4.10 Excel Output for Stresses ...................................................................................... 65 Figure 4.11 Comparison of 6’’ Span Results ........................................................................... 67 Figure 4.12 Comparison of 18’’ ............................................................................................... 67 Figure 4.13 Comparison of 30’’ ............................................................................................... 68 viii Figure Page Figure 4.14 Flow Chart of Computer Code for Improved Execution Speed ............................ 71 Figure 4.15 Linear Interpolation I ............................................................................................ 74 Figure 4.16 Flow Chart for Linear Interpolation I.................................................................... 76 Figure 4.17 Linear Interpolation II ........................................................................................... 78 Figure 4.18 Flow Chart for Linear Interpolation II .................................................................. 79 Figure 4.19 Excel Input for the Automated Program ............................................................... 81 Figure 4.20 Comparison of 6’’ Span Case with the Modified Code ........................................ 82 Figure 4.21 Comparison of 18’’ Span Case with the Modified Code ...................................... 82 Figure 4.22 Comparison of 30’’ Span Case with the Modified Code ...................................... 83 Figure 4.23 Slack Edge ............................................................................................................. 84 Figure 4.24 Slack Edge Criteria ............................................................................................... 85 Figure 4.25 Beams with Shear .................................................................................................. 86 Figure 4.26 Moment in a Beam ................................................................................................ 87 Figure 5.1 Tapered Roller Profile ............................................................................................. 92 Figure 5.2 Beisel Tapered Roller Model [12] .......................................................................... 93 Figure 5.3 The Figure for Tapered Roller ................................................................................ 96 Figure 5.4 Tapered Roller FE Wrinkle Model ......................................................................... 96 Figure 5.5 The Flow Chart for Tapered Roller ......................................................................... 99 Figure 5.6 Calculating Moment and Critical Taper ................................................................ 100 Figure 5.7 Excel Input of Tapered Roller ............................................................................... 102 Figure 5.8 Wrinkles Due to Taper, 92 ga Polyester, L=10” ................................................... 104 Figure 5.9 Wrinkles Due to Taper, 92 ga Polyester, L=20” ................................................... 104 Figure 5.10 Wrinkles Due to Taper, 92 ga Polyester, L=30” ................................................. 105 Figure 5.11 Wrinkles Due to Taper, 92 ga Polyester, L=40” ................................................. 105 Figure 5.12 Wrinkles Due to Taper, 56 ga Polyester, L=10” ................................................. 106 Figure 5.13 Wrinkles Due to Taper, 56 ga Polyester, L=20” ................................................. 107 Figure 5.14 Wrinkles Due to Taper, 56 ga Polyester, L=30” ................................................. 107 Figure 5.15 Wrinkles Due to Taper, 56 ga Polyester, L=40” ................................................. 108 Figure 6.1 Troughs and wrinkles due to a circular hole in the web [41] ................................ 110 Figure 6.2 Distance Between the Hole and the Tangent Line ................................................ 111 Figure 6.3 Mallya’s Model ..................................................................................................... 112 ix Figure Page Figure 6.4 The Figure for the Hole Model Traveling Between Rollers ................................. 113 Figure 6.5 The Web Hole FE Wrinkle Model ........................................................................ 114 Figure 6.6 Meshing the Hole Region...................................................................................... 115 Figure 6.7 Excel Output of Hole Region ................................................................................ 116 Figure 6.8 The Flowchart for Circular Discontinuity ............................................................. 118 Figure 6.9 The Section of Web over the roller ....................................................................... 120 Figure 6.10 Comparison of Results for Hole .......................................................................... 121 Figure 6.11 Excel output of square ......................................................................................... 122 Figure 6.12 Excel output of equilateral quadrangle ............................................................... 122 Figure 6.13 Comparison of Circle, Equilateral Quadrangle and Square Shaped Holes ......... 123 Figure 6.14 Flow Chart for Automating Hole Code ............................................................... 125 Figure 6.15 Excel Interface of Circular Void Excel VBA Code ............................................ 127 Figure 7. 1 Non uniform Web Models (COSMOS) ............................................................... 130 Figure 7. 2 Critical CMD Stresses Developed When the Non uniform Area is 3” Away From the Roller............................................................................................................... 131 Figure 7. 3 A Web with a Non uniform Circular Shaped Region Travelling Between Rollers132 Figure 7. 4 The Model for the Non uniform Hole Shaped Material ....................................... 133 Figure 7.5 Meshing the Non uniform Region......................................................................... 134 Figure 7. 6 Troughs and Wrinkles .......................................................................................... 135 Figure 7. 7 The Flowchart for Non uniform Webs ................................................................. 137 Figure 7. 8 Wrinkles Due to Non uniformities, 2r = 3”, t’ = 0.0008” .................................... 139 Figure 7. 9 Wrinkles Due to Non uniformities, 2r = 2”, t’ = 0.0008” .................................... 139 Figure 7. 10 Wrinkles Due to Non uniformities, 2r = 3”, t’ = 0.0001” .................................. 140 Figure 7. 11 Wrinkles Due to Non uniformities, 2r = 2”, t’ = 0.0001” .................................. 140 Figure 7. 12 Wrinkles Due to Non uniformities, L = 3”, t’ = 0.0005” ................................... 141 Figure 7. 13 Excel Interface of Non uniform Excel VBA Code ............................................ 143 Figure 8. 1 Meshing the Model .............................................................................................. 147 Figure 8. 2 The Flowchart for the Convergence Check ......................................................... 149 Figure 8. 3 Comparison of 6” Span Case after the Convergence Check ................................ 150 Figure 8. 4 Comparison of 18” Span Case after the Convergence Check .............................. 151 Figure 8. 5 Comparison of 30” Span Case after the Convergence Check .............................. 151 x Figure Page Figure 8. 6 Converging Result for a Specific Case for Misaligned Roller Case .................... 152 Figure 8. 7 Converging Result for a Long Narrow Web ........................................................ 153 Figure 8. 8 Converging Result for a Wide Web ..................................................................... 153 Figure 8. 9 Converging Result for 0.002” Thick 50” Wide Web ........................................... 154 Figure 8. 10 The Flowchart for the Tapered Roller ................................................................ 156 Figure 8. 11 Converging Results for 20” Case ....................................................................... 157 Figure 8. 11 Converging Results for 30” Case ....................................................................... 158 Figure 8. 11 Converging Results for 40” Case ....................................................................... 158 Figure 8. 14 Converging Result for a Specific Case for Tapered Roller Problem ................. 159 Figure 8.15 Converging Result for a 60” Long Web ............................................................. 160 Figure 8.16 Converging Result for a 80” Long Web ............................................................. 160 Figure 8. 17 Convergence Check for the Void Case .............................................................. 163 Figure 8. 18 Result of Modified Code for the Case at Figure 6.15 ........................................ 164 Figure 8.19 Convergence Result for a Wide Web for Circular Void Code............................ 165 Figure 8.20 Convergence Result for a Long Web for Circular Void Code ............................ 165 Figure 8. 21 Result of Modified Code for the Case at Figure 7.15 ........................................ 167 Figure 8. 22 Result of Modified Code for a Wide, Long Web ............................................... 168 1. CHAPTER I 1.1 INTRODUCTION A web is any flexible thin material. Webs are made in continuous production processes. The webs are often stored in wound roll form, since this is the only convenient means of storing long lengths of flexible material. A production roll of polyester film that is 48 gages (0.00048”) in thickness might be 4 feet in diameter. If the web is wound on a 6” diameter core there can be nearly 60 miles of web length in the roll. That polyester web could be nearly 300 inches wide when made, but costumer rolls to be shipped could be only 6 inches to 60 inches in width depending on the products to be made. Thus a particular web might be 60 miles long, 60 inches wide, and 0.00048 inches in thickness. Plastic films, papers, foils, and thin metal sheets are examples of webs. Web handling can be defined as all processes employed during the transportation of webs. Cutting, coating, slitting, printing, laminating and drying are some of the processes that add value to the web. During these processes, web materials must travel around several rollers in the process equipment. 2 During the transportation of webs trough process machinery, compressive stresses are induced in webs and these compressive stresses may cause instabilities. If these instabilities occur in free spans, which are the sections of web between two rollers, they are called troughs and if they occur at rollers, they are called wrinkles (Figure 1). Wrinkles and troughs result in loss of value and quality of the webs. For instance, if wrinkles may cause permanent damages to the webs the result is wasted web and highly priced downtime of web lines. The processes such as laminating, printing on paper, metalizing of films or coating require web to be planar. Troughs may affect these processes and like wrinkles the result is waste of material and time. As a result, wrinkles and troughs are two important engineering problems in the web handling industry. Figure 1. 1 Troughs and Wrinkles There are three types of wrinkles. Machine direction wrinkles, cross machine direction wrinkles and shear wrinkles. The machine direction [MD] is the direction of the web travel through the web process machine and the cross machine direction [CMD] is the direction that is perpendicular to the machine direction (i.e. across the web width). MD wrinkles occur 3 because of the compressive forces in the CMD. CMD wrinkles are produced during winding or unwinding processes due to interlayer slippage that may be the result of air entrainment [1]. Shear forces due to the imperfect rollers, misaligned rollers and non uniform webs may result in shear troughs and wrinkles. Shear wrinkles can result from misaligned or tapered rollers [1, 2]. Shear wrinkles can be classified as regime one wrinkles and regime two wrinkles. Roller Misalignment or Taper Web Tension Regime I Regime II Planar, Troughed, or Slack Edge Spans Wrinkles Figure 1. 2 Regime I and Regime II Wrinkles Regime one wrinkles are diagonal wrinkles which occur because of the presence of a lateral shear force in the web. The shear force might be due to a misaligned or tapered roller. The 4 shear forces in a web which are necessary to cause troughs or wrinkles are very small. The misalignments and tapers which cause the shear forces are unintentional. For regime one troughs and wrinkles to occur there must be sufficient traction between the web and the rollers to react the shear forces. The shear forces occur because when traction is adequate the web will attempt to gain normal entry to the downstream roller. The concept or law of normal entry of a web to a roller is attributed to Lorig [3]. Regime two wrinkles are dependent on traction and velocity between the web and the roller. In this case, if traction is insufficient, CMD compressive stresses cannot develop in the web. Frictional restraint forces between the web and the roller are required to sustain wrinkles. Increases in web tension result in increased normal forces between the web and roller and act to decrease the air films which develop between web and roller. Air films can develop between webs and rollers due to the hydrodynamic entrainment of air. Increased web tension and normal force act to decrease the thickness of the air films. The air films result in decrease in traction between the web and rollers. Thus an increase in web tension will increase the potential frictional forces that would be needed to sustain a wrinkle on a roller. The normal entry rule may be violated as a result [3]. Good et al [4] showed that regime two wrinkles could be predicted. In this proposal it will be assumed that sufficient traction exist between webs and rollers such that only regime one wrinkles need to be considered. Wrinkles are defined as buckling of webs around rollers. Webs have very small resistance to compression in free spans. Webs may withstand more compressive stress on the rollers than free span, because the web has greater stability in the form of a cylindrical shell than it has a flat plate. The critical CMD compressive stress needed to wrinkle the web can be predicted using classical cylindrical shell buckling expressions. 5 Previous research by Webb[9], Beisel[12] and Mallya[15] have shown that web wrinkles on rollers can be predicted by analyzing the compressive CMD stresses that form in the web as it contacts rollers. These compressive stresses arise due to troughs that have already formed in the free span prior to the rollers. The troughs themselves may be the result of many disturbances that exist in web lines. Misaligned rollers, tapered rollers and web non uniformities (such as holes) are examples of such disturbances. Thus, it could be stated that web wrinkles on rollers are the result of two instabilities. A disturbance is first required to induce a trough instability in the free span. The first instability has now occurred. The web that was planar in the free span has now troughed. After the troughs appear CMD compressive stresses will arise in the web entering the roller downstream of that span. When the troughs first appear, the CMD stress in the web entering the roller may be small. Thus, whatever disturbance produced the trough initially may have to become yet larger before sufficient CMD compressive stress can be generated in the web on the roller to produce an MD wrinkle. At this point, the second instability has occurred. The web on the roller that earlier had the shape of a sector of a cylindrical shell has now buckled. Thus, it has been proven that the prediction of wrinkles upon a roller involves: (1) An instability (troughs) occurring in the span upstream of the roller. (2) As whatever disturbance increases that produced the troughs in step (1), a post buckling analysis must be undertaken. It will bill seen that as the disturbance increases that a CMD compressive stress will arise in the web on the downstream roller. (3) As the disturbance yet increase further the CMD compressive stress in the web on the 6 downstream roller is also increasing. As this CMD stress increases it will finally surpass the cylindrical shell buckling stress. At this point, with sufficient friction between the web and the roller, a wrinkle will form in the web on the roller. The purpose of this study is to develop efficient computational tools that can accomplish the analyses required in step (1), (2), and (3). 7 2. CHAPTER II 2.1 LITERATURE REVIEW The literature has been reviewed and the findings will be broken into two sections. First the basic theories of membrane instability will be reviewed. Second, those studies which examine web instability will be reviewed. At the close of this chapter a final summary of the findings will be included and a statement of proposed research will be presented. 2.2 Theory of Membrane Instability Wagner [18] prepared a treatise on sheet metal girders with very thin webs. Probably this study is the earliest investigation of the mechanics of wrinkling membranes. He worked to develop the structural method of sheet metal girders. His methods were based on the 8 assumption of the low stiffness in bending of the metal web. He worked to explain the behavior of the thin metal webs in beams carrying a shear load well in excess of the initial buckling value. He proposed tension field theory. Tension field theory helps to analyze flexible structures that can support tension, but have no ability to resist compression. This theory was further developed by Reissner [23], SteinHedgepeth [19] and Mikolas [20, 21]. MillerHedgepeth [5] and Miller et al. [6, 22] adopted this theory to the finite element method. Tension field theory can be applied to the web lines due to the fact that web lines can support tension but cannot carry compression and also web materials are flexible structures. Stein and Hedgepeth [19] suggested a particularly useful approach concerning partly wrinkling membranes. This work is a seminal work in this field. They derived a theory to predict the stresses and deformations of stretched membrane structural components for loads under which part of the membrane wrinkles. Their theory was based on the basic assumption that a membrane has no bending stiffness and because of this can carry no compressive stress. They applied their theory to inplane bending of a stretched rectangular membrane, a pressurized cylinder, and to the rotation of a hub in a stretched circular membrane. They presented stresses and deformations in equation form for the wrinkled and unwrinkled regions. The membrane they considered is elastic, isotropic, has no bending stiffness, and cannot carry compressive stress. In their work they considered average strains and displacements of the wrinkling material rather than detailed deformations of each wrinkle. In terms of the wrinkling equations given, their theory was limited in the sense that average strains must be small compared with unity. They started to investigate the wrinkling region 9 by looking the principal stresses. They used a criterion that if both principal stresses are zero, the membrane is unloaded and thus will not wrinkle. Their criteria for a wrinkled membrane is one principal stress must be zero and other principal stress is non zero and tensile. The nonzero principal stress may be assumed to act along the crest of the wrinkle. The approach that was developed by Stein and Hedgepeth [19] was further developed and applied by Mikolas [20, 21]. He presented experiments and analysis for the wrinkling behavior of stretched membranes under the influence of a torque loading through an attached hub. He found that theory and experimental results were in a very good agreement. The work done by Stein and Hedgepeth [19] and Mikolas [20, 21] were closed form solutions and did not involve the finite element method. The principal stress criterion that they used is employed in our current approach. Miller and Hedgepeth [5] developed a new algorithm for finite element analysis based on the same assumptions and field equations after finding some critical disadvantages connected with the Stein and Hedgepeth [19] approach. This work may be the most important study in this field. In their algorithm the element stiffness is dependent on the current principal strains. Wrinkling membrane elements can have either taut behavior, wrinkled behavior or slack behavior. In taut behavior both principal stresses are larger than zero, in wrinkled behavior one of the principal stresses is greater than zero the other is equal to zero and in slack behavior both principal stresses are equal to zero. In their algorithm, in the first load step all elements are assumed to taut behavior. In the consecutive steps, element behaviors 10 are calculated with respect to strain states of previous steps. In other words, the decision on the stress state is made using the criteria based upon principal strains. In the algorithm they apply load step by step and they continue to solve for a particular load step until the convergence is achieved for that load step. In the proposed research, the algorithm that is employed is similar to Miller and Hedgepeth’s algorithm. Miller et al. [6, 22] investigated the algorithm further. They presented the efficiency and accuracy of the algorithm by applying it to the problems Stein and Hedgepeth [19] studied. They described the algorithm more detailed in these studies. They also described how Miller and Hedgepeth [5] derived the Dtaut, Dslack and Dwrinkled constitutive matrices in detail. 11 2.3 Applications of Instability Analyses to Webs Shelton [17] studied the steering effects of downstream misaligned roller. He modeled web span as a beam. His work helped us to justify boundary conditions which should be associated with a straight uniform web approaching a misaligned roller when there is sufficient friction to enforce the kinematic boundary conditions. Figure 2.1 Boundary Conditions for Misaligned Roller He determined that the concept of normal entry that had been used in the drive belt industry can also be used in web guiding systems. He stated that Vi and θi can be arbitrary set equal to zero without effecting the relative lateral deformation within the span. θj is the misalignment of the downstream roller at j, and due to the law of normal entry the slope of the web at j will be θj. Shelton proved that the final boundary condition was the downstream moment being zero. Thus Mj(L) = Vj’’(L) = 0. In our problem these boundary conditions will be used while modeling the instability of a web due to a downstream 12 misaligned roller. Gehlbach, et al. [32] proposed buckling criteria for troughed webs in a free span. They showed experimental verification for downstream misaligned roller case. In their work, isotropic web properties were considered. Gopal and Kedl [33] were the fist people who used finite element analysis and a commercial FE code to study trough formation in the web span between rollers. They modeled a web span by using triangular plate elements and they used ABAQUS FE commercial code. They were successful to predict deformations of a web due to the misaligned roller. Benson et al. [8] developed a finite element model for wrinkling analysis and called this code FEWA. They used this code to make calculations in their paper. Their aim was to better understand the conditions which cause wrinkle formation. They worked to predict locations where troughs would form and predict magnitudes of compressive stresses. They compared some of their results with the results from the nonlinear version of commercial code ABAQUS. In their code instead of using Wagner’s tension field theory they used the tension field theory that Wu [24, 25, and 26] introduced first. In this method, it is assumed that out of plane deflection relieves compressive stress across a wrinkle and that there is an associated strain with this deflection. Roisum [28, 29, 30, and 31] described wrinkling phenomena in detail. He explained 13 wrinkling, air entrainment, tension control, roller design, problems associated with profile variations, why and how wrinkles form, the types of wrinkles, troubleshooting techniques. He explained the importance of the problem due to the material cost and waste for the producers. While doing this he explained the subject in a simple way rather than an academic way. Shelton [14] worked on buckling of webs. He modeled a web as a buckled plate or a shell. He studied buckled wavelengths of webs in free spans and on rollers. He predicted the wavelength of the buckled form and he compared these results with experimental data. Good et al. [4] worked on velocity independent and velocity dependent wrinkles. Velocity dependence occurs when velocities are high enough and web tension is low enough to allow sufficient air entrainment between webs and rollers that Shelton’s boundary conditions [17] are no longer valid. They found that velocity dependent wrinkles can be avoided by using enough traction and suitable web line conditions. Hashimoto [34] worked on the studies done by Beisel [2, 10], which is a theoretical model for predicting the web wrinkling due to the misaligned roller. The theoretical model was established upon the experimental results. The experiments which he did were for non uniform webs with different Young’s modulus in MD and CMD directions. He compared these experimental results with the model and he verified the accuracy of the model. 14 Jones and McCann [35] studied wrinkling of webs on rollers and drums. They used a new variational analysis based on the method of RayleighRitz. They modeled the wrinkling as a continuous sine wave in a web on a roller or a drum. They reported that a shell buckling on a rigid support in an outward mode is similar to wrinkling on rollers. Beisel [2, 10] studied wrinkling phenomena due to the misaligned roller. In his study, he considered orthotropic material properties. He showed that while wrinkle formation would require post buckling analysis, prediction of trough formation can be expressed by closed form solution. Papandreadis [7] employed finite element methods to predict troughs in the webs. He studied effects of several parameters on the amount of lateral contraction of the web. He examined the effects of web material properties (Poisson’s ratio, modulus of elasticity) and web geometry (various lengthtowidth ratios), web thickness, loading conditions (tensile loading, combination of tension and shear forces) on the wrinkling phenomenon. He used a finite element code, named NASTRAN (Nasa Structural Analysis), to analyze the buckling of panels and the resulting shapes (like wavelength of the corrugations, number of waves). Webb [9] was the first person who tried to couple the behavior of the web in the free span to the web on the roller. Probably, the most important finding in his work was that CMD compressive stresses were forming in the web on the roller due to the troughs that had already formed in the free spans. He used quadrilateral elements within the commercial finite element code COSMOS to predict wrinkles due to the downstream misaligned roller. 15 He observed that there is a linear relationship between wrinkle θ and trough θ . He performed experiments to measure the deflection as a result of a misaligned roller. He found that the misalignment required to form a wrinkle was roughly twice that to form a trough over a wide range of parameters. He tested a wide range of web material properties, web thickness, and span length to width ratio (often called span ratio). He tried to find a relationship between the width reduction of the web and web buckling. He modeled the web crossing the roller using regular elastic elements. He used wrinkling membrane elements for the web in free span. In Fig.2.2 Webb’s approach to the problem is presented. Figure 2.2 Webb’s Model In Fig.2.3 the boundary conditions and the applied loads in Webb’s model is shown. 16 Force Prescribed Displacement Coupled DOFs Fixed DOF X Y Figure 2.3 Webb’s Model: Boundary Conditions and Applied Loads Timoshenko [11] showed the axial buckling stress for a sector of a cylinder is the same per unit circumferential length as that of a complete cylinder. He found that ycr σ is 3 (1 2 ) ycr h E R σ ν = − ⋅ ⋅ − (2.1) The web wrapped around the roller can be assumed to be sector of a cylindrical shell. Shelton [14] studied with the mechanics of buckling and he hypothesized that the web wrapped around the roller is shows similarity to an internally pressured thin cylinder vessel. He found that the tension in the web performs like the hoop stress in a pressurized cylinder and that a 17 pressure results between the web and roller. He discussed that Eq.2.1 may indeed be appropriate to use in modeling a web buckling on a roller. Good and Beisel [10] developed an instability criterion for orthotropic web in shell form in order to predict wrinkling of webs: 3(1 ) x y ycr x y t E E R v y σ = − − (2.2) where R is the radius of the cylinder, t is the thickness of the web, Ex and Ey are the elastic modulus in x and y directions respectively, and v is Poisson’s ratio. Beisel performed many tests of aluminum soda cans which were near perfect in geometry. He found that the buckling stress approached that given by Timoshenko’s expression and concluded that the earlier disagreement being due to shell imperfections must have been correct. Thus expressions 2.1 and 2.2 appear to be applicable to sectors of web transiting rollers. If the compressive stresses in the web on the roller reach the ycr σ , the web will buckle. Webb [9] increased the shear forces used to simulate the misalignment of the roller until the compressive stress induced in the web on the roller reaches the ycr σ value in (Eq.2.2). When the critical compressive CMD stress was reached, he determined the rotation of the downstream roller from the displacements output by the finite element code. The experimental critical misalignment of the roller and the rotation of the web at the entry of the 18 roller computed using finite element model did not match well. He proposed that the experimental critical misalignment of the roller cr θ is the sum of the critical angles predicted by the model for troughing, cr ,trough θ and for wrinkling, cr ,wrinkle θ in verifying his model. This was found later to be incorrect by Beisel [12]. Webb’s work was carried forward by Beisel [12]. Beisel studied troughing and wrinkling due to roller misalignment, roller taper, and roller crown. There is some similarity between the effects of roller misalignment and taper on troughs and wrinkles. In both cases the web is steered laterally in the machine and a shear stress results in the web. The misaligned roller induces the misalignment angle and a lateral deformation at the misaligned roller. The tapered roller induces a lateral deformation and a bending moment at the tapered roller. Beisel’s method of modeling the problem with COSMOS was different from the Webb’s method of modeling the problem. Beisel achieved good agreement with his models and the experimental results. His model and his results will be introduced in Chapter IV and Chapter V. Swift [39] examined steering of drive belts. He worked with the concept of a couple developed in a web approaching a tapered or crowned roller and the resulting steering of the web. He suggests the minimum amount of taper or crown which should be employed to control the web with minimum interference of stresses in the web. He gave experimental results to support his suggestions for corrective measures. 19 Shelton [44] discussed misaligned and tapered rollers, and he mentioned multiple web span interactions and moment transfer. In his thesis, Beisel [12] compared his model with Shelton’s model for the critical taper required to cause troughs to form in a web span approaching a tapered roller. Good and Beisel [45] worked on the formation of troughs due to tapered rollers. They attempted to determine whether the procedure that was employed for misaligned rollers would be applicable to the case of downstream tapered rollers. This work was a part of Beisel’s thesis [12]. Brown [46] presented a new method for modeling the elastic behavior of webs conveyed over rollers. He worked on lateral displacement of a misaligned roller and lateral displacement of a tapered roller. He suggested two modifications of the web boundary conditions. One was a generalization of the normal entry rule and the other was the addition of what he named the normal strain rule. With a numerical partial differential equation solver, he solved two dimensional plane stress equations and compared the results with earlier models. Shimizu et al. [47] and Shimizu [48] worked on plates which have holes and are subjected to tensile load. They used the finite element method in their work. They investigated the effects of aspect ratios (height/width) and shapes of holes to the k which is the buckling coefficient 20 of the plate described by Timoshenko and Gere [11]. They determined that with that curvatures on corners of holes have little effect in improving the tension buckling strength and that the buckling coefficient increases corresponding to the increasing aspect ratio. ElSawy and Martini [49] studied the effects of plate aspect ratio (height/width), hole size, hole location and loading ratio on the buckling coefficient k of rectangular perforated plates subjected to uniform end compression in the x direction and compression or tension in the y direction. Mallya [15] was the first person to examine the effects of holes in process machinery in web handling industry. He applied Beisel’s method of modeling to web wrinkles due to circular and elliptical discontinuities in the web. He studied the behavior of webs with voids traveling over a roller. He compared experimental results and finite element model results that he modeled using commercial FE code COSMOS. He studied elliptical voids and circular holes in terms of generating wrinkles. His FE model was similar in form to Beisel’s FE model with respect to boundary conditions and using five panels. His model and results will be introduced in Chapter VI. Kara [16] also used similar modeling method to predict the occurrence of wrinkles due to length variation across the web width. He attempted to find the critical conditions that would induce wrinkles. He heated the center of the web during his experiments to achieve length variation. He also used COSMOS FE commercial code to model his case, and used experimental findings to confirm his FE model. 21 Relevant books include, ‘Theory of Elastic Stability’ written by Timoshenko and Gere [11], ‘Introduction to Finite Elements in Engineering’ written by Chandrupatla and Belegundu [36], ‘Finite Elements and Approximation’ written by Zeinkiewicz and Morgan [37] and Visual basic Excel for Dummies [27]. Also ‘Finite Element Analysis Class Notes’ from Good [38] is used for finite element part of the study. 22 2.4 Summary The tension field theory can be applied to webs in web lines because the web can support tension but cannot carry compression. Wagner [18] proposed tension field theory. His theory was further developed by Raissner [23], Stein and Hedgepeth [19] and Mikolas [21, 22]. Miller and Hedgepeth [5] developed a new algorithm for finite element analysis based on this field equation. This algorithm is the most common algorithm used to examine wrinkles. Lorig [3] developed the normal entry rule for a web approaching a roller. Shelton [17] found that the moment in a web is zero when it approaches a misaligned roller. In the web handling area, Gopal and Kedl [33] were the first to employ a commercial finite element code to study trough formation in the web span between the rollers. Webb [9] attempted to model wrinkle formation due to a misaligned roller. He used partly wrinkling membrane elements while modeling the web span. He used the commercial finite element code called COSMOS while modeling his work, and shell buckling criteria to determine whether wrinkling occurred. Beisel [12] made the most recent attempts to model wrinkle formation due to the misaligned roller and the tapered roller. He modeled the web between the rollers by wrinkling membrane elements and linked it to the classical shell buckling criteria as Webb [9] did. He studied webs approaching misaligned rollers and tapered rollers. He used the commercial finite element code COSMOS to model these cases. He compared his model results with his experimental results .These results showed good agreement. Mallya [15] was the first person who investigated the effects of voids on the stability of 23 webs. He performed experiments and modeled circular and elliptical discontinuities in webs in web lines. He also used the commercial finite element code COSMOS in his models and his model results and test results matched well. 2.5 Statement of Proposed Research In the literature no citations were found for the wrinkling instability of moving webs that did not involve the use of commercial finite element codes for solution. Use of commercial finite element codes by novice users to solve nonlinear problems associated with web instability is difficult. The objective of the proposed research will be to develop user friendly finite element codes that will solve nonlinear instability problems associated with strain state dependent material properties and boundary conditions of moving webs. This code will be unique and will have economic value by helping minimize web material losses as described in the introduction. 24 3. CHAPTER III 3.1 FINITE ELEMENT EQUATIONS In this chapter, the finite element equations will be described in detail. Displacement, strain and stress equations, the element stiffness matrix, meshing, banded matrix and their relations will be studied. A discussion of the solution method for cases in which the elasticity matrix [D] is not constant will be given. 3.2 Two Dimensional Four Node Quadrilateral Elements In this study, two dimensional four node quadrilateral elements are used. In this section the equations and properties of four node quadrilateral elements will be given briefly. 25 The stiffness matrix for quadrilateral elements can be found from the strain energy term in a continuum. 1 2 T V U = ∫ σ ε dV (3.1) In finite elements we consider the total strain energy to be sum of the strain energies from each element. Eq.3.2 can be obtained by replacing dV by t dA in Eq.3.1, where t is the uniform thickness of an element. 1 2 T e e e U =Σt ∫ σ ε dA (3.2) The small strain displacement relations for two dimensional problems can be written as: x y xy u x v y u v y x ε ε ε γ ∂ ∂ ∂ = = ∂ ∂ ∂ + ∂ ∂ (3.3) Where u and v are the deformations in the x and y directions respectively. In twodimensional fields, the displacement components at each point in the domain of the finite element can be represented as functions of two coordinate directions. u = [u(x, y), v(x, y)]T (3.4) 26 For the general quadrilateral element shown in Fig.3.1, the nodal displacement vector is q = [q1,q2 ,q3 ,...,q8 ]T X Y u P(x,y) 4 2 1 3 q1 q2 q3 q4 q6 q8 q5 q7 v Figure 3.1 Four Node Quadrilateral Element The finite element method uses concept of shape functions to develop interpolations systematically. According to the concept, the shape functions must be developed for the master element. The master element is defined in natural coordinates ( ξ ,η ) and has a square shape in the natural coordinate system (Fig.3.2). The Lagrange shape functions are 1 2 3 4 N , N , N and N . Shape functions take the value of unity at the node where they are defined and the value zero at the 27 other nodes. For example 1 N takes the value one of unity at node one and takes the value zero at the other nodes. ξ η Figure 3. 2 Quadrilateral Elements Inξ ,η Space (Master Element) At the edges ξ = +1 andη = +1 1 N is equal to zero. So, 1 N must be a function like, 1 N = c (1−ξ ) (1−η ) (3.5) c is a constant that can be determined easily. Since 1 N is equal to one at node one, where ξ = −1 andη = −1.If we put these values at Eq.3.5, 1 = c (1− (−1)) (1− (−1)) (3.6) yields c =1/ 4 .Finally 1 N can be written as, 1 1 (1 ) (1 ) 4 N = −ξ −η (3.7) 28 By using the same procedure N1, N2 , N3 and N4 can be written as, 1 2 3 4 1 (1 ) (1 ) 4 1 (1 ) (1 ) 4 1 (1 ) (1 ) 4 1 (1 ) (1 ) 4 N N N N ξ η ξ η ξ η ξ η = − − = + − = + + = − + (3.8) Now, the shape functions can be used to interpolate the displacement at any point within the domain of the element using the equations: 1 1 2 3 3 5 4 7 1 2 2 4 3 6 4 8 u N q N q N q N q v N q N q N q N q = + + + = + + + (3.9) if N is, 1 2 3 4 1 2 3 4 0 0 0 0 0 0 0 0 N N N N N N N N (3.10) then the displacement can be written as: [ ]{ } 1 2 3 1 2 3 4 4 1 2 3 4 5 6 7 8 0 0 0 0 [u] 0 0 0 0 q q q u N N N N q N q v N N N N q q q q = = = (3.11) with the help of shape functions a point in the element can be described as, 29 1 1 2 2 3 3 4 4 1 1 2 2 3 3 4 4 x N x N x N x N x y N y N y N y N y = + + + = + + + (3.12) The shape functions can also be used to generate a map between the cartesian coordinates (x,y) and the natural coordinates (ξ ,η ). Since the same shape functions have been used to interpolate deformation within an element and to generate the coordinate map equations this is called an isoparametric formulation. The strain relationships (3.3) require derivatives with respect to cartesian coordinates. We currently have the deformations u and v defined as shape functions, which are functions of the natural coordinates ξ andη , multiplied by nodal deformations that are constants. So, to determine the strains in cartesian coordinates we must first relate the derivatives of deformations in natural coordinates (ξ ,η ) to derivatives in cartesian coordinates(x,y). If a displacement function in x , y coordinates is u = u(x , y) then this function can be considered to be an implicit function of ξ andη asu = u[x(ξ ,η ) , y(ξ ,η ) ] . Differentiation due to the chain rule, u u x u y x y u u x u y x y ξ ξ ξ η η η ∂ ∂ ∂ ∂ ∂ = + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = + ∂ ∂ ∂ ∂ ∂ (3.13) If we define the Jacobian matrix as, J x y x y ξ ξ η η ∂ ∂ ∂ ∂ = ∂ ∂ ∂ ∂ (3.14) We can rewrite Eq.3.13 as, 30 =J u x y u u x x u x y u u y y ξ ξ ξ η η η ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ (3.15) It seems that the Jacobian can transform derivatives in cartesian coordinates to derivatives with respect to natural coordinates. By using Eq.3.8, Eq.3.12 can be written as, 1 2 3 4 1 2 3 4 1 1 1 1 (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) 4 4 4 4 1 1 1 1 (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) 4 4 4 4 x x x x x y y y y y ξ η ξ η ξ η ξ η ξ η ξ η ξ η ξ η = − − + + − + + + + − + = − − + + − + + + + − + (3.16) With the help of Eq.3.16 Jacobian term Eq.3.14 can be written as, 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) J 4 (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) x x x x y y y y x x x x y y y y η η η η η η η η ξ ξ ξ ξ ξ ξ ξ ξ − − + − + + − + − − + − + + − + = − − − + + + − − − − − + + + + − (3.17) this equation can be written as, 11 12 21 22 J J J J J = (3.18) Now the Jacobian can be inverted and rewrite Eq.3.15 to produce derivatives that are related to strains: J 1 u u x u u y ξ η − ∂ ∂ ∂ ∂ = ∂ ∂ ∂ ∂ (3.19) 31 by using Eq.3.17, Eq.3.19 can be written as 22 12 21 11 1 det J u u x J J u J J u y ξ η ∂ ∂ ∂ − ∂ = ∂ − ∂ ∂ ∂ (3.20) by following the same procedure v displacements can be obtained as, 22 12 21 11 1 det J v v x J J v J J v y ξ η ∂ ∂ ∂ − ∂ = ∂ − ∂ ∂ ∂ (3.21) The equation ∂x ∂y = det J ∂ξ ∂η has a proof in reference [36]. By using Eq.3.3, Eq.3.20, Eq.3.21 and defining A matrice as, 22 12 21 11 22 11 22 12 0 0 1 0 0 det J J J A J J J J J J − = − − − (3.22) The strain displacement relations can be written as 22 12 21 11 22 11 22 12 0 0 1 0 0 det J x y xy u u u x u u J J v J J A y v v J J J J u v y x v v ξ ξ ε η η ε ε γ ξ ξ η η ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ − ∂ ∂ ∂ = = = − = ∂ ∂ ∂ − − ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ ∂ ∂ (3.23) with the help of 3.17 G matrice can be defined as, 32 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 1 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 4 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) G η η η η ξ ξ ξ ξ η η η η ξ ξ ξ ξ − − − + − + − − − + + − = − − − + − + − − − + + − (3.24) and using the displacement vector Eq.3.9 the derivatives of u and v can be written in the natural coordinates as a function of the nodal deformations, (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 1 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 4 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) u u q G q v v ξ η η η η η ξ ξ ξ ξ η η η η ξ ξ ξ ξ ξ η ∂ ∂ ∂ − − − + − + ∂ − − − + + − = = ∂ − − − + − + ∂ − − − + + − ∂ ∂ (3.25) If B is defined as B=AG, by using Eq.3.23 and Eq.3.25 strain term can be written as, 1 2 3 22 12 4 21 11 5 22 11 22 12 6 7 8 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 0 0 1 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 0 0 4det J 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) q q q J J q J J q J J J J q q q η η η η ξ ξ ξ ξ ε η η η η ξ ξ ξ ξ − − − + − + − − − − + + − = − − − − + − + − − − − − + + − B q = (3.26) two dimensional constitutive relations will be used to relate stress to strain σ =Dε and now, σ =D Bq (3.27) 33 By using Eq.3.26 and Eq.3.27 the strain energy term 1 2 T e e e U =Σt ∫ σ ε dA , becomes, 1 1 1 1 1 [ det ] 2 1 2 T T e e T e e U q t B DB J d d q q k q ξ η − − = = Σ ∫ ∫ Σ (3.28) where k is 8x8 element stiffness matrix: 1 1 1 1 e T det e k t B DB J dξ dη − − = ∫ ∫ (3.29) These integrals can be evaluated by using numerical integration methods. The Gaussian approach will be considered for this purpose. 3.3 Numerical Integration by Gaussian Approach By integrating in natural coordinates the bounds of integration are much simplified. In cartesian coordinates the ybounds will be functions of x and the xbounds will be functions of y. In natural coordinates our bounds are from 1 to1 for both ξ andη . Series can be used to take the integrals (Eq.3.29).By using the Gaussian quadrature approach; integration can be evaluated using weights and points. These points are also called Gauss points. 1 1 1 2 2 1 ( ) ( ) ( ) ( ) n n I f ξ dξ ω f ξ ω f ξ ω f ξ − = ∫ ≈ + +L + (3.30) 34 Eq.3.30 is an example of n point approximation. Here ω1,ω2K ωn are weights and 1 2 , n ξ ξ K ξ are the Gauss points. If a one point formula is employed, the integral becomes: 1 1 1 1 I f (ξ ) dξ ω f (ξ ) − = ∫ ≈ (3.31) If a two point formula is employed, the integral becomes: 1 1 1 2 2 1 I f (ξ ) dξ ω f (ξ ) ω f (ξ ) − = ∫ ≈ + (3.32) The finite element method naturally incorporates some error as a numerical approximation. The complex continuums were modeled with many finite elements with simple shape functions to represent the element deformations (Eq.3.8). Hence, it would be undesirable to incorporate additional error in the integration of stiffness terms. It is desirable that the integration be exact. In a one point formula two parameters are considered 1 1 (ω and ξ ) . Suppose that our integration is required be exact when f (ξ ) a polynomial of is order one. So, suppose f (ξ ) is a function 0 1 f (ξ ) = a + a ξ . If we select 1 1 ω = 2and ξ = 0 Gaussian quadrature will yield an exact result. In a two point formula there are four parameters to choose 1 2 1 2 (ω ,ω ,ξ and ξ ) .Suppose that f (ξ ) must be exact for a cubic polynomial, 2 3 0 1 2 3 f (ξ ) = a + a ξ + a ξ + a ξ . The error term will be, 1 2 3 0 1 2 3 1 1 2 2 1 Error (a a ξ a ξ a ξ ) dξ [ω f (ξ ) ω f (ξ )] 0 − = ∫ + + + − + = (3.33) 35 Solution yields four nonlinear equations and they have the unique solution, 1 2 1 2 1 1/ 3 ω ω ξ ξ = = = − = − (3.34) As a result by using two Gauss points and by using the values in Eq.3.34 a cubic expression or less can be integrated exactly. By increasing the number of Gauss points, different weights (ω) and different locations (ξ ) can be found. In the FE algorithm that is developed in this study two Gauss points are used. Two Dimensional Integrals The equation 3.29 for our stiffness terms involves twodimensional integrals. So we need to extend Gaussian quadrature to the two dimensional integral form: 1 1 1 1 I f (ξ ,η) dξ dη − − = ∫ ∫ (3.35) If I is in a form like Eq.3.35, I can be written as, 1 1 1 1 1 1 1 [ ( , ) ] [ ( , ) ] ( , ) n i i i n n j i i j j i n n i j i j i j I w f d w w f w w f ξ η η ξ η ξ η − = = = = = ≈ ≈ ≈ ∫ Σ Σ Σ ΣΣ (3.36) Stiffness matrix (Eq.3.29) is twodimensional integral. The product of BTDBdet J is quadratic in terms of ξ andη .So two point Gauss Quadrature yield an exact result. It is 8x8 36 matrix so it has 64 elements. Each term must be calculated by using Eq.3.36. If ( , ) ( T det ) e ij f ξ η = t B DB J putting two for n in Eq.3.36 yields, 2 2 1 1 1 1 2 1 2 2 1 2 1 2 2 2 ( , ) ( , ) ( , ) ( , ) ij k ≈ w f ξ η + w w f ξ η + w w f ξ η + w f ξ η (3.37) where 1 2 1 1 2 2 1 1 1 3 3 w w ξ η and ξ η = = = = − = = (3.38) After input of these weights and Gauss points into the Eq.3.37 kij can be found as, 1 1 1 1 1 1 1 1 ( , ) ( , ) ( , ) ( , ) 3 3 3 3 3 3 3 3 ij k ≈ f − − + f − + f − + f (3.39) ξ η 1 1 ( , ) 3 3 1 1 ( , ) 3 3 1 1 − ( , ) 3 3 − − 1 1 ( , ) 3 3 − 1 2 3 4 Figure 3.3 Two Points Gauss Quadrature using 2x2 rule. 37 3.4 The Global Stiffness Matrix and Matrix Banding The element stiffness matrix is an 8x8 matrix and it has 64 elements. After calculating all stiffness matrices for all elements, the global stiffness matrix need to be formed from the element stiffness matrices. In Fig. 3.4 the plate is divided into n elements, this procedure is called meshing. For every element each of the stiffness terms ( , ) ( T det ) e ij f ξ η = t B DB J is evaluated using Eq.3.29. While calculating the terms of global stiffness matrix, if a node is only used by one element, its stiffness terms should be placed in the global stiffness matrix directly (like Fig.3.4 node 1). If a node is used by two elements (like Fig.3.4 nodes 3,5..) , stiffness terms for this node from two elements must be added to form the stiffness terms of the global stiffness matrix for that node. Lastly if a node is used by four elements (like Fig.3.4 nodes 4,6…) stiffness terms for that node from four elements must be added to form the global stiffness matrix. Figure 3.4 A Simple Model for Global Stiffness Matrix 38 Figure 3.5 A Quadrilateral Element and Nodal Displacements In Eq.3.40 and Eq.3.41 the stiffness matrices for element one and element two displayed in Fig.3.4 are shown 1 1 1 1 1 1 1 1 11 12 13 14 15 16 17 18 1 1 1 1 1 1 1 1 22 23 24 25 26 27 28 2 1 1 1 1 1 1 33 34 35 36 37 38 3 1 1 1 1 1 44 45 46 47 48 4 1 1 1 1 5 55 56 57 58 1 1 1 6 66 67 68 1 1 7 77 78 1 8 88 q q q q q . q q q k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k sym k k k k k k 1 2 3 4 5 6 7 8 F F F F F F F F = (3.40) 39 2 2 2 2 2 2 2 2 11 12 13 14 15 16 17 18 5 2 2 2 2 2 2 2 22 23 24 25 26 27 28 6 2 2 2 2 2 2 33 34 35 36 37 38 7 2 2 2 2 2 44 45 46 47 48 8 2 2 2 2 9 55 56 57 58 2 2 2 10 66 67 68 2 2 11 77 78 2 12 88 q q q q q . q q q k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k sym k k k k k k 5 6 7 8 9 10 11 12 F F F F F F F F = (3.41) After global stiffness matrix is formed, all displacements are calculated by using Gaussian elimination method. After finding displacements, strains and stresses for all elements can be calculated. In the algorithm that is developed the user decides the mesh density. 1000 or more elements may be needed for long spans simulations. The resulting stiffness matrix will be large and much computational time will be required for solutions. Computational times can be large because a multistep solution will be required where loads are slowly increased and the [D] matrices updated after each load step. Thus the size of the stiffness matrix becomes important because the system of updated equations will be solved many times. The stiffness matrix in our problem is a symmetric matrix. Instead of using the whole stiffness matrix the banded form of the stiffness matrix can be employed and reduce the computation time. To explain the form of the banded matrix, assume that the plate in Fig.3.4 is meshed with only two elements. In Eq.3.30 and 3.31 the stiffness matrixes of two elements were given. So, the global stiffness matrix for the plate will be: 40 1 1 1 1 1 1 1 1 11 12 13 14 15 16 17 18 0 0 0 0 1 1 1 1 1 1 1 22 23 24 25 26 27 28 0 0 0 0 1 1 1 1 1 1 33 34 35 36 37 38 0 0 0 0 1 1 1 1 1 44 45 46 47 48 0 0 0 0 1 2 1 2 1 2 1 2 2 2 2 2 55 11 56 12 57 13 58 14 15 16 17 18 1 2 1 2 1 2 66 22 67 23 68 24 k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k + + + + + + + q1 q2 q3 q4 q5 2 2 2 2 25 26 27 28 q6 1 2 1 2 2 2 2 2 q7 . 77 33 78 34 35 36 37 38 1 2 2 2 2 2 q8 88 44 45 46 47 48 2 2 2 2 q9 55 56 57 58 q10 2 2 2 66 67 68 q11 2 2 77 78 q12 2 88 k k k sym k k k k k k k k k k k k k k k k k k k k k k k k + + + F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 = (3.42) The banded form of 3.42 is shown in 3.43. The bandwidth of the matrix has been reduced to 8. 1 1 1 1 1 1 1 1 11 12 13 14 15 16 17 18 1 1 1 1 1 1 1 22 23 24 25 26 27 28 0 1 1 1 1 1 1 33 34 35 36 37 38 0 0 1 1 1 1 1 44 45 46 47 48 0 0 0 1 2 1 2 1 2 1 2 2 2 2 2 55 11 56 12 57 13 58 14 15 16 17 18 1 2 1 2 1 2 2 2 2 66 22 67 23 68 24 25 26 2 k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k + + + + + + + 2 7 28 1 2 1 2 2 2 2 2 77 33 78 34 35 36 37 38 1 2 2 2 2 2 88 44 45 46 47 48 2 2 2 2 55 56 57 58 2 2 2 66 67 68 2 2 77 78 2 88 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 k k k k k k k k k k k k k k k k k k k k k k k k k + + + (3.43) Banded matrix equation solvers exist which helps to greatly reduce the computational time associated with solving set the set of equations during each load step. 41 To conclude, the element stiffness terms are calculated from Eq.3.29 where ( , ) ( T det ) e ij f ξ η = t B DB J . After developing element stiffness matrices, the global stiffness matrix is formed. After the total potential energy is formed, which is the sum of the strain energy terms from Eq.3.28 and the work potential of internal and external forces, the theory of minimum total potential can be used. The resulting systems of equations that must be solved have the form: K q = F (3.44) This set of equations is arranged in the banded form discussed and then solved for the nodal deformations {q}. After finding displacements strains can be calculated from, ε = Bq (3.45) The stress and elastic strain components are related by a set of coefficients known as Generalized Hooke’s Law. This law can be written as, σ = Dε = D B q (3.46) Where D is equivalent elasticity matrix. 42 3.5 A [D] Matrix for Plane StressIsotropic Material and Partly Wrinkled Membrane Elements The D matrix is used while forming element stiffness matrixes and finding element stresses. The algorithm that is developed is nonlinear since the membrane D matrix is dependent on principal strains. In our models, elastic elements are used for the web on the rollers and for the exiting free span. Partly wrinkled membrane elements that were developed by Miller and Hedgepeth [5] are used for the web in test free span. The D matrix can be explicitly stated for any material. For the plane stress state for an isotropic material, stress strain relations can be written as: 1 ( ) 1 ( ) 2 (1 ) x x y y y x xy xy v E v E v E ε σ σ ε σ σ γ τ = − = − = + (3.47) By solving for the stresses in Eq.3.47 the D matrix can be found: 2 1 0 1 0 1 (1 ) 0 0 2 v E D v v v = − − (3.48) In taut membrane behavior, both principal stresses are greater than zero and the web cannot trough or buckles out of plane. In taut behavior, the two inplane principal strains may both be greater than zero or 1 ε can be greater than zero and 2 ε less than zero as long as the ratio 43 (ε 2 /ε1 ) is greater than the negative value of Poisson’s ratio ( −v ) . The [D] matrix given in Eq.3.48 is used to relate stress to strain in (1) membrane elements exhibiting taut behavior and (2) in linear elastic elements, used to model the web upon rollers. For the linear elastic elements there are no conditions placed upon the use of this [D] matrix, the principal stresses and strains can take on positive and negative values. When unsupported membrane elements have a first principal strain 1 ε less than zero this infers that the second principal strain is also negative and less than 1 ε from the rules we use to order principal strains. In this case the membrane elements exhibit a slack or unstressed behavior between stress and strain and the [D] matrix is null. In wrinkled membrane behavior the first principal stress 1 σ is greater than zero and the second principal stress 2 σ is zero. In terms of strain this behavior is represented with a first principal strain 1 ε greater than zero but now the second principal strain 2 ε is always negative. Not only 2 ε is negative, but the ratio of ( 2 1 ε /ε < −v ). For a simple membrane in uniform tension with no lateral constraint we would expect a lateral contraction governed by the expression 2 1 ε = −vε . We would expect this lateral contraction occur while the web remained planar, no buckling would be expected. However if 2 ε becomes more negative than 1 −vε we would expect a compressive 2 σ stress to have developed but since membranes by definition can withstand no compressive 2 σ stress without buckling, we would assume this element has entered the wrinkled state. For the wrinkled membrane positive nonzero principal stress can be supposed to act along the wrinkle. In a wrinkled membrane element in longitudinal tension, if 1 σ is the nonzero positive principal stress, the longitudinal direction 44 will be observed to be parallel to the wrinkles. By the property of Mohr’s circle it is known that 1 2 2 2 x y R O ε ε ε ε − = + = (3.49) where R is the radius of Mohr’s circle for strains and O is the distance from the origin of the coordinate system to the center of the circle. Consider a thin flat membrane in a state of plane stress in an xy coordinate system. Let the principal axes I and II be rotated through an angle α relative to the xy axes (Figure 3.6). 2α xy γ x ε y ε ε Figure 3.6 Mohr’s Circle for Plane Strain From the Fig.3.6 45 1 2 1 2 cos 2 2 cos 2 2 x y R R ε ε ε α ε ε ε α + = + + = − (3.50) Substituting the value of R (Eq.3.49) into Eq.3.50 yields 1 2 1 2 1 2 1 2 ( ) cos 2 2 2 ( ) cos 2 2 2 x y ε ε ε ε ε α ε ε ε ε ε α − + = + + − = − (3.51) Eq. 3.51 yields 1 2 1 2 (1 cos 2 ) (1 cos 2 ) 2 2 (1 cos 2 ) (1 cos 2 ) 2 2 x y ε ε ε α α ε ε ε α α = + + − = − + + (3.52) If P = cos 2α , then Eq.3.52 becomes 1 2 1 2 (1 ) (1 ) 2 2 (1 ) (1 ) 2 2 x y P P P P ε ε ε ε ε ε = + + − = − + + (3.53) Also xy γ can be calculated from Fig.3.6. 1 2 1 2 2 sin 2 2( ) sin 2 ( ) sin 2 2 xy xy γ R α ε ε γ α ε ε α = − = = − (3.54) If Q = sin 2α then Eq.3.54 becomes 1 2 ( ) xy γ = ε −ε Q (3.55) 46 As a result, the usual strain transformation equations yield, 1 2 1 2 1 2 (1 ) (1 ) 2 2 (1 ) (1 ) 2 2 x y xy P P P P Q Q ε ε ε ε ε ε γ ε ε = + + − = − + + = − (3.56) where 1 2 1 2 cos 2 sin 2 x y xy P Q ε ε α ε ε γ α ε ε − = = − = = − (3.57) Within a wrinkled region usual elasticity equations do not apply. Instead, the assumption of negligible bending stress in the membrane yields the stress state as: σ1 = Eε1 ; σ 2 = 0 (3.58) Or 1 1 1 1 1 1 1 1 1 1 (1 ) (1 ) 2 2 (1 ) (1 ) 2 2 1 1 2 2 x x y y xy xy P P E P P E Q QE σ σ σ ε ε ε σ σ σ ε τ σ τ ε = + = + = − ⇒ = − = = (3.59) Expressing stresses in terms of strains in the form of σ = Dε is desirable for numerical analysis. D is 3×3 matrix, ( , , ) T ( , , ) T x y xy x y xy σ = σ σ τ and ε = ε ε γ . Because the problem is statically determinant within the wrinkled region (for example 2 σ = 0 independent of the 47 valuesε1 and ε 2 ) D is singular and many possible representations for D are possible. If λ plays role of Poisson’s ratio then 2 1 λ = −ε /ε . In a wrinkled region λ is always greater Poisson’s ratio [19]. In some points λ can take the value 1. So, choosing D matrix similar to Eq.3.48 is not a right choice because of the 1/(1−λ 2 ) term. Consider D matrix is like a b c b d e c e f (3.60) If Eq.3.60 substitute into σ = Dε and impose Eq.3.56 and Eq.3.59 yields (1 ) (1 ) 2 0 0 0 (1 ) (1 ) (1 ) 2 0 0 0 0 0 (1 ) 0 (1 ) 2 0 (1 ) 0 (1 ) 0 (1 ) 2 0 0 0 0 (1 ) 0 (1 ) 2 0 0 (1 ) 0 (1 ) 2 0 P P Q a P P P Q b P P Q c P E P P Q d P P Q e Q P P Q f + − + − + − + − − = − + − + − − + − (3.61) Solutions for the elements of D matrix are not unique because the coefficient matrix in Eq.3.61 is not singular. If b=0 selected and replaced into Eq.3.61 it is found that (1 ) ; (1 ) 2 2 ; 4 4 E E a P d P E E c e Q f = + = − = = = (3.62) The resulting D matrix is the matrix which is presented by Miller and Hedgepeth [5] for partly wrinkled membranes. D matrix is: 48 2(1 ) 0 0 2(1 ) 4 0 1 P Q E D P Q Q + = − (3.63) where P and Q stated at Eq.3.57. As a result, in our models elastic elements and wrinkled membrane elements are used. For the elements in the web region on the rollers, the D matrix (Eq.3.48) is used to relate stress to strain. The elements in the web span have strain dependent D matrices as explained below. Although the slack behavior is possible in some applications it is not possible in a nonlinear formulation which employs an incremental force method. Once edge slackness begins in a web using a force method, convergence to an equilibrium solution is not possible. Convergence would be possible using an incremental displacement solution. Since this research is focused on applications where edge slackness does not occur, an incremental force solution was acceptable. 1 σ 1 σ 1 σ 1 σ 2 σ 2 σ 2 σ = 0 1 σ = 0 2 σ = 0 1 1 2 0 <ε and υε > −ε 1 1 2 0 <ε and υε < −ε Figure 3.7 Behaviors of Wrinkling Membrane Elements [12] D matrices are defined for all allowable behaviors (taut, wrinkled and slack). In Fig 3.7 behaviors of wrinkling membrane elements can be seen. A useful algorithm for choosing D 49 matrix may be expressed as [5]: 1 1 1 2 ; 0 ; 0 ; s w T D D D D and D D otherwise ε ε υε ε = < = > < − = (3.64) Where D matrices are defined as, 0 s D = (3.65) For slack behavior, as 2 1 0 1 0 1 0 0 (1 ) / 2 T E D υ υ υ υ = − − (3.66) For taut behavior, and as 2(1 ) 0 0 2(1 ) 4 1 W P Q E D P Q Q Q + = − (3.67) For wrinkled behavior where P = cos 2α ,Q = sin 2α and as stated in Eq.3.57. These algorithms will be used to establish the code that will be explained in the following chapters. 50 4. CHAPTER IV 4.1 MISALIGNED ROLLERS In this chapter the work done of Beisel [12] will be briefly reviewed. The algorithm developed will be explained step by step. Also the code that implements the algorithm will be described. Then the measures taken to decrease CPU time and to automate the code will be discussed. Finally a new slack edge criterion for misaligned rollers will be developed. 4.2 Beisel’s Method for Modeling Wrinkles Due to the Misaligned Rollers Beisel developed a method to model wrinkle formation due to a downstream misaligned roller. Like Webb he used commercial FE code COSMOS. He used wrinkling membrane elements in free spans to allow troughs to form. In Fig.4.1 Beisel’s wrinkling model for 51 misaligned roller is shown. Figure 4.1 Beisel’s FE Wrinkle Model for Misaligned Roller Contrary to Webb, he used five panels of elements rather than three panels. The first panel represents the web on the upstream roller. The second panel represents the free span where wrinkling membrane elements are used. The third panel represents the web on the downstream misaligned roller. He employed the fourth and fifth panels to enforce desired boundary conditions. First, he applied tension to the web when he reached the desired tension load he began to apply shear force to the model as shown in Fig. 4.1. Beisel employed this five panel model for the following reasons: A. The asymmetric shear forces allowed him to model the zero moment boundary condition at the misaligned roller found by Shelton [17]. B. The fourth panel was modeled using regular elastic elements. This was done because the fourth panel acts to increase the bending stiffness of the elastic elements in panel three. The failure sequence of events was: 1. Troughs form at a critical angle of misalignment given by Beisel’s previous 52 work [2]. In the finite element code the troughs are modeled by elements that may assume the wrinkled membrane behavior described in Chapter III. This happens whenever σ 2 < 0 . In a real web this does not occur until the 2 σ stress becomes more negative than the buckling stress in the free span ycr σ ( 3(1 2 ) x ycr h E a π σ σ ν = − − ; h is thickness, and a is the distance between two rollers [12]). The onset of troughs can be predicted with linear buckling theory with closed form expressions developed by Beisel [12]. To predict wrinkling of the web on rollers requires nonlinear analysis since the entering free span has already buckled in the form of troughs. 2. After troughs form CMD compressive stresses begin to appear in the elastic elements in panel three that border panel two. As the shear forces and the associated misalignment is increased further, these CMD compressive stresses become more negative and finally approach the value in Eq. 2.1, at which point wrinkles are eminent. The elastic elements in the panel four restrict the bending in panel three due to the troughs that have formed in panel two. Beisel increased the shear forces until the critical compressive stress given by Eq. 2.1 was induced in the linear elastic elements at the entry of the misaligned downstream roller. Then, he concluded that the rotation of the nodes at the entry of the downstream roller should be equal to the angle misalignment in the roller. Beisel and Webb ran similar experiments to determine the onset of wrinkle formation due to the misaligned downstream roller. Beisel compared his results with these experimental results. He modeled a polyester web with a thickness of 0.00092 in (92 gages). The web 53 parameters for this web were a Young’ Modulus of 712000 psi, a Poisson’s ratio of 0.3, a web width of 6’’ and again the thickness was 0.00092’’. The rollers had a radius of 1.45’’ and ycr σ from the Eq.2.1 was about 270 psi. He modeled 6’’; 18’’ and 30’’ span lengths and compared his results with the experimental findings. In Fig.4.2Fig.4.4 the comparison of model and experimental data is presented. 0 0.002 0.004 0.006 0.008 0.01 0.012 0 5 10 15 20 25 30 35 40 45 Tension (lbs) θcr (rad) . . experimental data FE wrinkle model Figure 4.2 L = 6’’ Span Results 0 0.005 0.01 0.015 0.02 0.025 0 5 10 15 20 25 30 35 40 45 Tension (lbs) cr (rad) . experimental data FE wrinkle model Figure 4.3 L = 18’’ Span Results 54 0 0.01 0.02 0.03 0.04 0 5 10 15 20 25 30 35 40 45 Tension (lbs) θcr (rad) . experimental data FE wrinkle model Figure 4.4 L = 30’’ Span Results As seen from these three charts, the results from FE model and experimental results show good agreement. So Beisel’s model was successful in estimating wrinkle formation due to misaligned roller for a typical web span. In his study he compared his results with Webb’s results. He claims that his model yields better results than Webb’s model. He also claimed that the assumption that Webb proposed ( cr θ is the sum of cr,trough θ and cr,wrinkle θ ) is not true. Beisel achieved good agreement with experimental results without relying upon Webb’s assumption. 4.3 A New Algorithm for Predicting Wrinkles Due to the Misaligned Rollers As mentioned in the introduction, there are three types of wrinkles: MD wrinkles, CMD wrinkles and shear wrinkles. Shear wrinkles can occur due to roller imperfections such as 55 misaligned rollers or tapered rollers. The goal of this new algorithm is to codify and automate the analysis that Beisel perfected using a commercial finite element code (COSMOS). The term “shear wrinkle” resulted from the realization that these troughs and finally wrinkles were the result of shear forces in the web. Both the misaligned roller and the tapered roller induce shear in the web. Beisel [12] was successful to develop a method for predicting web wrinkles on rollers by using membrane elements described by Miller and Hedgepeth [5]. He applied this method to the prediction of wrinkles due to misalignment in rollers, tapered rollers and crowned rollers and he confirmed his results with laboratory tests. Beisel used commercial finite element code COSMOS to apply this method. In this method while the elements representing the web on rollers are modeled with elastic elements, the web in the free spans are modeled with wrinkle membrane elements. These elements cannot react compressive stresses and they can be in one of three states. These states include taut web, wrinkled web and slack web. In the taut web state, the elements can resist tensile stresses in both principal directions. In the wrinkled web state, membrane elements can withstand tensile principal stress in one direction and zero stress in the other principal direction. In the slack web state, the elements can carry no stresses in any direction. In this algorithm, forces are increased in to the model step by step. In the first step, all elements are modeled with elastic elements. After the first step, the principal strains for each element in free spans are calculated and stored. By using principal strains, code will select which of D matrices given below for three states will be used for the next load step. 56 1 1 1 2 ; 0 ; 0 ; s w T D D D D and D D otherwise ε ε υε ε = < = > < − = (4.1) where D matrices are defined at Eq.3.653.67. After convergence is obtained in each step, the compressive CMD stresses in the linear elastic elements are reviewed. If those stresses remain greater than the Timoshenko shell buckling stress (Eq.2.1) the shear force would be increased. If the CMD stresses in these elements became more negative than the Timoshenko shell buckling stress (Eq.2.1), the shear force would need to decrease and a bracketing method would be employed to determine what shear force would produce a negative a negative CMD stress essentially equal to the Timoshenko buckling stress. Once accomplished, the misalignment or taper that induced that level of shear force would be determined. By using the method explained above, a finite element code will be developed in Excel VBA (Visual Basic Excel). This code can be executed in any PC with Excel installed without need of a commercial FE code license and it allows users to analyze the misaligned roller case using a simple Excel based interface. The advantage of this code will be that users will not need any linear or nonlinear finite element background to execute the code. They do not have to know the kinetic and kinematic boundary conditions for misaligned rollers. The inputs will include parameters such as web tension, web width, span length, roller diameter, Poisson’s ratio and elastic modulus. When executed the code will automatically form a finite element mesh based upon the inputs with elastic quadrilateral elements representing the web supported by rollers and with wrinkle membrane quadrilateral elements representing the web 57 in the free span. The first code will implement boundary conditions for a web approaching a misaligned roller. Other boundary conditions will be studied later. Beisel [2, 12] and Webb [9] studied the misaligned roller case. The boundary conditions that they used are first proposed by Shelton [14] and then Good et al. [4]. They considered the web span as a beam. A classic beam is one which the web span length would be ten times longer than the width. Shear effects become “important” when L 10 w < . Tension becomes “important” when the lateral deformations become large. “important” in this context means that these effects have sizable influence on the lateral deformations of the web. The boundary conditions that are used by Beisel [2, 12], Webb [9], and others will be used in this model. The validity of using these boundary conditions was verified by comparison to experimental results by these authors. The normal entry condition of a web approaching a roller was enforced using coupling equations which enforce multipoint constraints. Lines of adjacent nodes crossing a roller in panel one and panel five had their CMD displacements coupled. Each adjacent line of nodes was coupled separately and in this way Poisson contraction of the web could occur unimpeded. There was no coupling of nodes in the web in contact with the misaligned roller. Since the moment in the web in the vicinity of the misaligned roller is small or zero the deformations of nodes are nearly that associated with a rigid body rotation. This results in the normal entry condition being satisfied in the web at the misaligned roller without enforcing a 58 multipoint constraint. The lines of nodes in the CMD at the exit of panel one and entry of panel five were each coupled in machine direction displacement. This procedure was done to ensure the maximum moment in the free web spans occurred at the border between span one and span two and the border between span four and span five. The system that is modeled is shown in Fig .4.5. Figure 4.5 The System That is Modeled The system of five panels is shown in Fig.4.5 modeled in Fig.4.6. The coupling discussed earlier is also shown. 59 Figure 4.6 Misaligned Roller FE Wrinkle Model The model is divided into five sections. The first panel represents the web on the upstream roller. The coupled nodes in this panel are used to enforce normal entry and exit on upstream roller. The second panel represents the entry web span to the misaligned roller. Here wrinkling membrane elements are used to simulate web behavior which allow troughs. Different from Beisel’s model at the first attempt the fourth panel is also modeled with wrinkling membrane elements. Shear forces are applied to the web on the upstream and downstream rollers to simulate the shear, moment, and lateral deformations of a web passing over a misaligned roller. The third panel represents the web on the downstream misaligned roller. A central node is fixed in the MD and CMD directions to prevent rigid body translations of the model. Rigid body rotation is prevented by the coupling of the CMD deformations of the lines of nodes crossing panel one and panel five. The fourth panel and the fifth panel elements and boundary conditions help ensure the zero moment boundary condition at the misaligned roller. The flow chart for the program is shown in Fig. 4.7. 60 61 Figure 4.7 Flow Chart for New Algorithm 62 The stiffness terms for the elements on the upstream roller are calculated by using DT. Then these stiffness terms are assembled into the global stiffness matrix. The same procedure is followed for the other web regions on the rollers. For the web spans D matrix for the elements will be formed differently. At the first load level and first iteration DT is used for all elements. Then, for all remaining load levels and for all iterations, the program selects one of the three D matrices from (Eq. 3.64) Dtaut, Dwrinkled or Dslack by evaluating the principal strains calculated in a previous step as explained in Eq. 3.6567. The nonlinearity for this case is due to the variable D matrices for the span elements as the shear loads increase. After selecting D matrices, elemental stiffness matrices are formed and the stiffness terms will be assembled into the global stiffness matrix using the same procedure as the web on the rollers. After the global stiffness matrix is formed, lines of nodes in the MD on the first and fifth panels are coupled in the y direction and the point at the center of the model is fixed in x and y directions as shown in figure 4.6. Next, the shear and traction forces are applied to the system. From the set of equations KQ=F the displacements {Q} can now be calculated. Strains {ε } and stresses {σ } can then be determined using the displacements. The strains are calculated usingε = BQ , stresses are calculated with the aid of σ = DBQ for all elements. After calculating strains, principal strains are also determined from the cartesian strains so that the proper D matrices can be selected for the next iteration or for the next load level. The shear and tension forces are applied incrementally in 20 steps. Five iterations are performed for each load step to allow P and Q to converge. As mentioned before if the 63 principal axes one and two are rotated through an angle α relative to the xy axes, P is the cosine of that angle and Q is the sine of that angle. For each new load step, the code first uses the principal strains calculated in the 5th iteration of the previous load step to determine the state of the wrinkling membrane for the 1st iteration. Then, the strains calculated in the first iteration for that load level are used to select D matrices for the second iteration. The procedure will continue in the same way. At each iteration, the strains calculated for the previous iteration are used to select D matrices till the maximum iteration number (5) is reached for that load step. The analysis proceeds in load steps with iterative analysis steps occurring within each load step. The iterative analysis steps are necessary to allow the values of P and Q to converge in the elements with the wrinkled behavior prior to moving to the next load step. In Fig.4.8 it shows how P and Q (Eq.3.57) behave with iteration for an element in Panel 2. Figure 4.8 Convergences of P and Q 64 Figure 4.9 Excel Input for the Program In Fig. 4.9, the screen where the web parameters are input to the code is shown. As seen from the Excel input screen, the user enters the parameters such as width of the web, a roll dimension (quarter circumference of the roller), a span length dimension, the thickness of the web, the elastic modulus of the web, Poisson’s Ratio of the web, and the web tension in units of traction (stress) in the x direction. Also the shear in y direction in units of traction (stress) to the Excel sheet is input. In the “m height elements” cell, user enters the mesh density along x direction. In “n1 roll elements” cell, user inputs the mesh density of the rollers along y direction. Similarly in n2 span elements cell user decides about the mesh density of the span elements along y direction. This is not the final form of input. In the next chapters the efforts for the final form of input will be discussed. 65 The output of the code developed to predict wrinkling due to downstream misaligned roller are σ x and σ y stresses for all elements on upstream roller, the downstream roller and in the web span between the upstream and the downstream roller are shown in Figure 4.10. Figure 4.10 Excel Output for Stresses For this case critical buckling stress by Eq.2.1 is about 270 psi. The marked row of stresses are the stresses of the elements at the entry of the downstream roller. The elements at the entry of the roller buckle first. After the output is displayed the minimum (most negative) stress in these elements should be compared with the critical shell buckling stress. If unequal, then the shear force should be increased or decreased until the minimum compressive stress 66 in the elements at the entry of the downstream roller reaches the critical shell buckling stress. In Fig.4.10, the maximum compressive stress (269.6) reached the Timoshenko shell buckling stress (270). At this instant, the angle displayed as an output in Figure 4.9. is the critical angle of misalignment of the downstream misaligned roller for the onset of wrinkling (0.0073 rad.). 4.4 Comparison of Results with Previous Works The code that was developed was executed for some cases defined by Beisel and the code results are compared with Beisel’s experimental and FE results. As mentioned in the previous chapter he used a commercial FE code called COSMOS and modeled the misaligned roller case. He performed his tests using 92 gage polyester web to verify his model. For this web, Young’s Modulus is 712000 psi, Poisson’s ratio is 0.3, the width of web is 6’’ and the thickness is 0.00092’’. The roller has a radius of 1.45’’. Eq.2.1 yields about 270 psi for the critical shell buckling stress for this case. Beisel was comparing the stresses at the nodes at the entry of the roller with the critical compressive stress predicted by Timoshenko shell buckling criteria. In this study, compressive stresses in the elements at the entry of the roller are used in the comparison instead of nodal stresses. Stresses at four Gauss points of the elements are computed and their averages are found as elemental stresses. The following graphs in Fig. 4.1113 show the experimental critical angle of misalignment at 67 the downstream roller, the angle predicted by Beisel’s model and by Yurtcu’s model for 6’’, 18” and 30’’ web span lengths. 0 0.005 0.01 0.015 0.02 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester,L=6" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu's Model(VBA) Figure 4.11 Comparison of 6’’ Span Results 0 0.005 0.01 0.015 0.02 0.025 0.03 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester, L=18" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu's Model (VBA) Figure 4.12 Comparison of 18’’ Span Results 68 0 0.01 0.02 0.03 0.04 0.05 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester,L=30" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu's Model (VBA) Figure 4.13 Comparison of 30’’ Span Results From the figures, it can be seen that the results of the code developed agree well with the lab test data and the results of the commercial code. 4.5 Improvement of Execution Time It was found that execution the code for long spans that twenty load steps with five internal iterations in each load step required a large amount of CPU time. For example, 2000 elements may be employed in a long span case. For every element we have an (8*8) stiffness matrix. Forming the global stiffness matrix for 2000 elements and solving it 100 times (five iterations in 20 load steps) required extensive CPU time. Although we use banded matrix in 69 our code long spans could require three hours to produce a result on a single core computer. Due to these long execution times the next focus was decreasing the CPU time. First the number of iterations was decreased from five to two, then from two to one, and the result was not significantly changed. Then step by step the load step increments were decreased from twenty to four, and still reasonable results were achieved. However, it was found that decreasing the load step increments to less than four caused the results to err dramatically. This problem is appears to converge with four load increments and one iteration within each load increment. It was mentioned in the previous section that Beisel modeled the second span with elastic elements, but that wrinkling membrane elements were used in this code. This means that we were calculating principal strains to select which D matrix would be used for every element in span two. Beisel ran taut elements in the second free span because the model was a 900 wrap case. For the misaligned roller case this subjects the upstream span to shear and the downstream span to twist. A web span will absorb a large amount of twist without forming negative 2 σ stresses. So Beisel assumed the D in all elements in the downstream span would remain taut. We forced D matrices to be DT in the second free span. This also provided better agreement between our code and other results. Changing element types from wrinkling membrane elements to elastic elements also helped to decrease CPU time. Our code produced reasonable results, and did this within seconds. After making these changes, the flowchart of our program will be like Fig.4.14. The name “WRINKLINGsystem” will be applied to the part of the new chart that begins after “Mesh Model with Quadrilateral Elements” and continues to “Load Level<4”. This name will be used while attempting to automate the code in Part 4.6. 70 71 Figure 4.14 Flow Chart of Computer Code for Improved Execution Speed 72 4.6 Automating the Code From Fig 4.9 it can be seen that the user enters parameters such as web width, roller wrap arc length, span length, web thickness, web elastic modulus, Poisson’s ratio, the number of elements across the web width, the number of elements down the span length, the number of elements across the roller, the web tension and shear. The research objective was to limit the inputs that could be comprehended by users with no knowledge of the finite element method. These would include only web width, roller wrap arc length, span length, web thickness, web elastic modulus, Poisson’s ratio and web tension. One of the inputs required by the code is web tension. Web tension is one of a very few parameters that are controllable in a web process machine. Thus it would be optimal for a chart to be produced for the user that shows how much misalignment is allowable as a function of web tension, rather than computing what misalignment in a roller is acceptable at one tension. The user can then decide to solve an instability problem by better aligning the rollers or by changing the web tension. Thus the user should enter only certain parameters and the code will determine the mesh parameters and shear force required to induce wrinkling. To automate the mesh parameters the code was executed several times to explore what mesh density was required to produce a converged result. Very reasonable results were obtained by dividing the width and span length per piece per dimension and dividing the rollers with the integer part of six times of onefourth of roller circumference. The mesh density is delimited for very long span lengths, very short span lengths, very long wide webs and very short wide webs. The meshing procedure and convergence check will be addressed more detailed at 73 Chapter VIII. The next step was determining how to automate the search for the shear force that was required to induce the wrinkle instability. The user inputs the web tension. A linear interpolation scheme was used to determine the shear force. One level of shear will produce a certain level of compressive y σ stress at the misaligned roller, which can be determined by the code. A second level of shear will produce another level of compressive y σ stress. Interpolation can be used to estimate the level of shear that will produce the shell buckling stress. It is an estimate because this is nonlinear analysis. That estimate can then be input to the code to help refine the actual shear level that will induce wrinkling. A slack edge criteria was used as a starting point. If a slack edge forms during the computations the code will fail. This is because increased shear will not result in increased compressive stress in the web on the roller, it will only increase slackness. For more information about slack edge criteria earlier work done by Good [40] can be visited. Web spans can be modeled as beams [12]. From EulerBernoulli beam bending theory the stress in the crosssection is: F My A I σ = ± (4.2) If this stress is equal to zero, a slack edge occurs. If our web has a thickness of t, width of w, span length of L and the applied traction in x direction is Tx the shear stress (Txy) for the web to be slack can be calculated from Eq.4.2 as: 0 2 x w M T I σ = = − (4.3) Where moment M is calculated from: 74 M = Txy ×w×t × L (4.4) If this value put into Eq.4.2 and do the calculations Txy (Tslack) will be found as 6 x xy T w T L × = (4.5) This was used as the starting value of Traction XY (shear) in the code (Fig.4.9).It found that for all cases computed the compressive stress in the first row of elements in the roller two was lower (less negative) than the value that was calculated from Timoshenko buckling criteria (Eq.2.1). Thus it found that if the shear value in Eq.4.5 was used as traction xy (Tslack) and half of it ( Tslack/2) one could be sure that there were two data points which would produce compressive stresses less than that is needed to buckle the web (Tcritical) (Figure 4.15). critical σ 2 σ 1 σ Figure 4.15 Linear Interpolation I 75 After finding two points we can estimate Tcritical from the tangent tanα : tan 2 1 2 ( / 2) critical Slack Slack Critical Slack T T T T σ σ σ σ α − − = = − − (4.6) from here Tcritical can be found as: 2 2 1 ( ) ( ( / 2)) critical Slack Slack ( / 2) Critical Slack T T T T σ σ σ σ − × − = + − (4.7) The algorithm for Linear Interpretation I can be seen in Figure 4.16. The code takes the value of web tension (traction x) from Excel input page and calculates F1, F2 and F3 values with the help of Tslack (Eq.4.5). Here F1 is equal to Tslack /2, F2 is equal to Tslack and F3 is equal to 2* Tslack .Than code runs for F1 and F2 values within the WRINKLINGsystem and calculates Sigma1 and Sigma2 values. By using Eq.4.7 code calculates Tcritical . The code checks whether Tcritical is larger than two times of Tslack . This was done because during the runs it was observed that if this was not done, the value of Tcritical increases dramatically because of the angle between two points. Taking this step helps control Tcritical value. Then the code runs WRINKLINGsystem for Tcritical value and finds SigmaS. The code continues this process until SigmaS is bigger than the Timoshenko buckling criteria (Sigmacritical). At the end of Linear Interpolation I, there is a F1 value which is less than the traction that is needed to buckle the web, and there is a F2 value which is more than the traction that is needed to buckle the web (Figure 4.17). 76 Figure 4.16 Flow Chart for Linear Interpolation I 77 The code for linear interpolation one will be: TRACTIONXY = (TRACTIONX * H) / (6 * (L1 + L2 + L1 + L2 + L1)) F2 = TRACTIONXY F3 = 2 * TRACTIONXY Do Until sigmaS > syrc FORCEXY = F2 * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigma2 = critical F1 = F2 / 2 FORCEXY = F1 * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigma1 = critical FS = ((sycr  sigma2) * (F2  F1)) / (sigma2  sigma1) + F2 If FS > F3 And control = 0 Then FS = F3 control = control + 1 End If FORCEXY = FS * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigmaS = critical F1 = F2 F2 = FS Loop 78 After establishing values F1 and F2 we can start Linear Interpolation II process. The new problem will look like Figure 4.17. critical σ 2 σ 1 σ Figure 4.17 Linear Interpolation II From Fig.4.17 Tcritical can be found as: 1 2 1 ( ) ( 2 1) critical 1 critical F F T F σ σ σ σ − × − = + − (4.8) Here the aim is to approach Tcritical value by changing the values of F1 and F2. The flowchart for Linear Interpolation II can be seen at Figure 4.18. 79 Figure 4.18 Flow Chart for Linear Interpolation II From Linear Interpolation I the values F1 and F2 are known. By using Eq.4.8 Tcritical was calculated and SigmaS was found. If SigmaS is greater than Sigmacritical we replace Tcritical with F2. If not we replace Tcritical with F1. This process was continued until SigmaS is 80 between the limits. In case of misaligned roller case the lower limit was set to 0.99*Sigmacritical and the upper limit was set to 1.01*Sigmacritical. The code for Linear Interpolation II will then be: Downlimit = 0.99 * sycr Uplimit = 1.01 * sycr Do Until Downlimit <= sigmaS And sigmaS <= Uplimit FORCEXY = F2 * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigma2 = critical FORCEXY = F1 * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigma1 = critical FS = ((sycr  sigma1) * (F2  F1)) / (sigma2  sigma1) + F1 FORCEXY = FS * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigmaS = critical If sigmaS > sycr Then F2 = FS ElseIf sigmaS < sycr Then F1 = FS End If Loop Linear Interpolation I, Linear interpolation II and WRINKLINGsystem can be seen in the 81 Appendix. After automating, Excel input of the code will look like Figure 4.19. Figure 4.19 Excel Input for the Automated Program Here the user is supposed to enter web width, span length, thickness, elastic modulus, Poisson’s ratio, roller radius and the web tension. Other parameters are calculated automatically. Meshing elements are calculated with Excel equations and shear traction xy required to induce wrinkles is calculated within the code as mentioned above. After execution the code provides the following output: total time of execution (67 seconds), the maximum compressive stress (266.8 psi) in the first row of elements and the roller misalignment angle (in degrees) that produced that compressive stress. For this case the sigma critical value calculated from Timoshenko buckling criteria (Eq.2.1) is 270 psi. So the interpolation scheme discussed has produced a compressive stress in the first row of elements 82 that is very close to the critical value. The results from the automated code are shown in Figures 4.2022 and the results are compared with Beisel’s commercial FE results and his test data. 0 0.005 0.01 0.015 0.02 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester,L=6" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu`s Modified Code (VBA) Figure 4.20 Comparison of 6’’ Span Case with the Modified Code 0 0.005 0.01 0.015 0.02 0.025 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester, L=18" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu`s Modified Code (VBA) Figure 4.21 Comparison of 18’’ Span Case with the Modified Code 83 0 0.01 0.02 0.03 0.04 0.05 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester,L=30" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu`s Modified Code (VBA) Figure 4.22 Comparison of 30’’ Span Case with the Modified Code The execution time of the modified code is much less than the previous version. In the previous version execution time was around two to three hours and code was able to give one result for the specified shear force. For the modified code the execution time is around five to ten minutes and the modified code finds the right shear force that will buckle the web by itself. Over the parameter ranges of these examples the modified code appears to mesh the problem adequately and yields good results at all tension levels. 84 4.7 A New Slack Edge Criteria for a Misaligned Roller In Figure 4.23, buckling region of a 92 gage polyester web can be seen. For this web Young’s Modulus is 725000 psi, Poisson’s ratio is 0.3, the width of web is 6’’, length of the web is 40” and the thickness is 0.00092’’. If 6000 psi MD stress applied to the system we are able to see wrinkles after the trough formation. If 2000 psi applied to the system after the formation of troughs slack edge occurs and we are not able to see wrinkles. Figure 4.23 Slack Edge The code was delivered to the Web Handling Research Center sponsors. One of the sponsors used the code with a low span (L/W) ratio. In that case W (width of the web) was five times 85 larger than L (length of the web). With that specific web tension that was applied to the web a slack edge was supposed to occur. But Excel VBA code was came up with a result which means wrinkles was occurred after the trough formation instead of a slack edge. The case was also modeled with COSMOS and the result was very similar with Excel VBA code. In the derivation of Eq.4.5 shear deformation was neglected. For L/W values of 0.2 it was obvious that shear deformation is important. Thus it was determined that a new slack edge criteria was needed that did account for shear deformation. Figure 4.24 Slack Edge Criteria The previous slack edge expression was derived based on Euler beam theory: ( ) x SlackEdge T a b Ehb θ = (4.9) 86 Here, a is length of the web, b is width of the web (Fig.4.23), h is thickness and E represents elastic modulus. To derive a new criteria the Timoshenko beam theory was used that incorporates both shear and tension stiffening. Figure 4.25 Beams with Shear If w is deformation, the slope of the beam isα , and can be found from: ' w wb ws w x x x α θ γ ∂ ∂ ∂ = = = + = + ∂ ∂ ∂ (4.10) Here θ is related to moment and γ is related to shear. ( ) ( ) s M x F x and x EI GA θ γ ∂ = = ∂ (4.11) The moment can be calculated from Eq.4.10 by using Eq.4.11 ( ) 1 " , , . x x s M x F Moment w EI GA x θ γ ∂ = = + = + ∂ (4.12) 87 Figure 4.26 Moment in a Beam At the boundary x = 0 , 0 0 x x S w F x GA γ = = ∂ = = ∂ (4.13) The moment M(x) = F(L − x) and F(x) = F , and substituting into Eq.4.12 yields ( ) " 0 F L x w EI − = + (4.14) Integrating yields 2 1 ' ( ) 2 F x w Lx C EI = − + (4.15) From Eq.4.13 1 x 0 S F C GA γ = = = so Eq.4.15 becomes 2 ' ( ) 2 S F x F w Lx EI GA = − + (4.16) Finally integration of Eq.4.16 yields, 88 3 2 2 ( ) 6 S F x Fx w Lx C EI GA = − + + (4.17) It is known that 0 0 x w = = . So 2 C becomes 0. We seek critical θ at x=L .From Eq.4.16, 2 1 ' ( ) x L critical 2 S L w F EI GA θ = = = + (4.18) Solving for the steering force F yields: 2 1 ( ) 2 critical S F L EI GA θ = + (4.19) This allows us to calculate the moment at x=0, 2 . 1 ( ) 2 critical S M F L L L EI GA θ = = + (4.20) The beam bending stress from moment can now be calculated as, max 2 1 1 2 ( ) 2 B critical S M w y L I L I EI GA θ σ = = + (4.21) The normal stress due to web tension can be calculated from, T T wb σ = (4.22) In taut behavior σ T >σ B if σ T =σ B the slack case occurs. So if σ T =σ B , 89 2 1 2 ( ) 2 critical S L w T L I wb EI GA θ = + (4.23) From here critical θ can be calculated as, 2 2 2 ( ) critical S T L I Lw b E GA θ = + (4.24) For a rectangular cross section 3 5 , 12 2(1 ) 6 S bw E I G and A bw v = = = + . Substituting into Eq.4.24 yields, 2 2 2 (1 (1 )) ( ) 5 critical TL w v w Ebw L θ = + + (4.25) In the calculations beginning from Eq4.10, L for span length, w for span width and t for thickness were used. If we replace them with a for length, b for width and h for thickness, it will be found, 2 2 (1 ( ) (1 )) ( ) 5 critical Ta b v b Ehb a θ = + + (4.26) If this result is compared with Eq. 4.9 it will be found that 2 2 (( ) (1 )) ( ) 5 Ta b v b Ehb a + (4.27) is the effect of shear force with Timoshenko beam theory. If v ≈ 0.3 then the term is 90 2 1 (1 ) 5 2 + v ≈ . So Eq.4.26 can be simplified as, 2 1 (1 ( ) ) ( ) 2 critical Ta b b Ehb a θ = + (4.28) If W >L (b>a) then shear effects are becoming significant. For instance if W=5L as it was in the case that we came up with earlier then, 1 (1 25) ( ) 2 13 ( ) criticalTimoshenko criticalEuler Ta b Ehb Ta b Ehb θ θ + ≅ ≅ (4.29) Eq. 4.26 supersedes Eq.4.9 and can be used for both long and short spans. 91 5. CHAPTER V 5.1 TAPERED ROLLERS In this chapter, a brief definition of tapered rollers, previous work modeling with tapered roller cases with commercial finite element codes will be reviewed. Then a new model based on using Excel VBA will be given. 5.2 Description of Tapered Rollers and Beisel’s Method for Modeling Wrinkles Due to the Tapered Rollers A tapered roller is defined as a roller with a linearly varying radius across its width [12]. Tapered rollers are commonly seen in the web handling industry. The process of roller manufacture will almost certainly result in rollers with a slight taper. These tapered rollers 92 are produced unintentionally; the use of tapered rollers in web lines may result in web damage. To solve this problem, machining techniques that involve feedback can be used, but this will be costly. Therefore, knowing the amount of taper that will not result in harm to the web process is beneficial for the industry. In Figure 5.1 a web approaching the tapered roller is shown. In the figure, the taper is somewhat overemphasized so that it can be seen. y r Roller radius Position across the web Average radius of roller Ro Figure 5.1 Tapered Roller Profile The radius of the roller at any point across the width is: r( y) = my + Ro (5.1) Here m is the slope of the roller. The velocity across the roller width is: ( ) ( ) ( ) o V y = r y ω = my + R ω (5.2) Here ω is the angular velocity of the roller. The average web velocity can be found by using average roller radius as: 93 Vavg = Roω (5.3) Variation in the velocity across the web width will cause a stress and strain upon the web: ( ) ( ) avg ( ) ( ) md md avg o o V y V my Emy y and y E y V R R ε σ ε − = = = = (5.4) These equations assume the web and roller achieve the same velocity at contact. The variation in stress across the web width causes a steering moment on the web: ( ) 2 2 2 3 2 2 12 b b j b b o o Emhy mEhb M y hydy dy R R σ − − − − = ∫ − = ∫ = (5.5) where m is roller taper, R0 is roller radius, E is Young’s Modulus, h is thickness and b is web width. Beisel [12] determined that wrinkle formation due to a downstream tapered roller is similar to the wrinkle formation due to a misaligned roller. As in the misaligned roller case, he used Timoshenko buckling criteria (Eq.2.1) to decide whether the web on the tapered roller will wrinkle or not. His wrinkle model for the tapered roller is shown in the Figure 5.2. Figure 5.2 Beisel Tapered Roller Model [12] Web span being studied Wrinkling membrane elements Nodes coupled along dotted lines, to provide normal entry, exit and travel across rollers Horizontal dotted lines coupled in y direction Vertical dotted lines coupled in x direction Uniform shear force (2fyj) applied at both ends of middle roller and end of web spans Web line tension applied at both ends of web Nodes coupled in y direction to 94 As in the misaligned case, he divided the model into five panels. The first panel is the web on the upstream roller. The second panel represents the web span being studied. The third panel represents the tapered roller. The fourth panel represents the web between the tapered roller and third roller, and the fifth panel represents the third roller. He locked the nodes in the first panel along horizontal lines. Each node on the line must move with the same y displacement as the rest of the nodes on the same line. This allows for Poisson contraction due to web line tension but requires a zero slope at the beginning of the troughed web span. He also locked the nodes at the exit of the upstream roller in the x direction along a vertical line to simulate the roller gripping the web in a noslip condition. He used wrinkling membrane elements in the second panel. Along the right edge of the second panel he locked the nodes along horizontal lines in the y direction for a very short distance. He did this to ensure the normal entry of the web to the tapered roller. He modeled the right hand side of the model to enforce the boundary conditions. He executed the model by first applying web line tension and then increasing shear force until the compressive stress across the first row of elastic nodes on the elastic wrinkling membrane element boundary reached the critical value predicted by Eq. 2.1. Then he calculated the moment associated with the last row of span elements by using Eq. 5.5. He calculated the critical taper that would induce wrinkles. He obtained the experimental results for the onset of wrinkles due to a downstream tapered roller. He compared his experimental results and model results for two materials. In the following chapters, we will compare our model results with his experimental results and with his model results. 95 5.3 A New Algorithm for Predicting Wrinkles Due to the Tapered Rollers As mentioned earlier, shear wrinkles can occur due to tapered rollers. Beisel [12] successfully modeled the tapered roller case by using the commercial finite element code COSMOS. As in the misaligned case he used wrinkling membrane elements while modeling the elements between the upstream roller and the tapered roller. The goal of our new algorithm related to the tapered roller is to codify and automate the analysis that Beisel did with the commercial finite element code COSMOS. Similar to the misaligned roller case a finite element code that calculates critical taper will be developed in Excel VBA by using wrinkling membrane elements. The boundary conditions will be similar to Beisel’s model boundary conditions. The problem is modeled material on the upstream roller, the upstream span, the web on the tapered roller, the downstream span, and finally material on the downstream roller. The system that is modeled is shown in Fig.5.3. 96 Figure 5.3 The Figure for Tapered Roller The system with a tapered roller shown in Figure 5.3 is modeled like the Figure 5.4 shown below. Figure 5.4 Tapered Roller FE Wrinkle Model 97 The new model is similar to the misaligned roller case. The model was divided into five panels. The first panel was the panel on the upstream roller. The second panel was representing the web between upstream roller and tapered roller. Here wrinkling membrane elements were used; in the rest of the model elastic elements were used. The third panel represents the tapered roller. The fourth and fifth panels help to achieve the no moment boundary condition at the tapered roller. Boundary conditions and loads were enforced which Beisel found to be appropriate for the tapered roller. The center is pinned to prevent rigid body motions. Multipoint constraints were applied on each row of nodes on the entering and exiting rollers. Along the right edge of the second panel (span one), Beisel locked the nodes along horizontal lines in the y direction for a very short distance to ensure normal entry. In contrast to his method, high mesh was used at the last row of the span one element (ten elements were used in the last row elements) and six points along the horizontal lines were locked in the y direction. This was decided after trying many ways to achieve normal entry to the tapered roller. The flowchart for the program is shown in Fig.5.5. It is similar to the misaligned roller flowchart so it will not explain it here in detail. 98 99 Figure 5.5 The Flow Chart for Tapered Roller 100 In the tapered roller flowchart the forces are applied in five time steps. After five steps from the first row of elements the moment is calculated and by using Eq. 5.5 critical taper is also calculated. Assume that the stresses of the first row of elements at tapered roller look like Fig.5.6. Here, one half of the width is demonstrated. hi Ai n ı i Figure 5.6 Calculating Moment and Critical Taper Moment is M = F h and / x σ = F A. From these two equations the moment of an element can be calculated from 1 n i i i i M σ A h = =Σ . If all moments are added the total moment of the first row of the elements of the tapered roller is found. If this value put into Eq .5.5 the critical taper for that specific case is found. At the end of the code, the code calculates total moment and taper for that case. The name “SystemWrinkler” will be applied to the part of the flow chart (Fig.5.5) that begins after “Mesh Model with Quadrilateral Elements” and continues to “Load Level <5” .This term will be used in the following chapter. 101 5.4 Automating the Code Similar to the misaligned roller case, it was aimed for the user to enter only certain parameters such as web width, roller wrap dimension, span length, web thickness, elastic modulus, Poisson’s ratio, and web tension, and run the code. The attempts to automate the code began with trying to automate the mesh. After running many cases, it was determined that reasonable results were achieved by using three elements per dimension for the web width, one element for every dimension of the web spans and the integer part of four times of the roll dimension for rollers. Excel equations were used to set these values. The meshing procedure and convergence check will be addressed more detailed at Chapter VIII. The flowcharts for Linear Interpolations I and II are used for automating the tapered roller case. Similar methods to the misaligned case were applied. Tslack and Tslack/2 (Eq.4.5) were used as a starting point for Linear Interpolation I. In Fig.4.16 and Fig.4.18, if we use the term SystemWrinkler instead of WRINKLINGsystem, the way the new flow chart works can be explained. After automating our tapered roller code the Excel input of the code will look like Figure 5.7. 102 Figure 5.7 Excel Input of Tapered Roller Here the user is supposed to enter web width, span dimension, thickness, elastic modulus, Poisson’s ratio, roller radius and the web tension. After entering these parameters with the help of the Excel equations the mesh parameters, one fourth of the circumferences of the rollers (L1 roll dimension in Fig.5.7) and the stress calculated from Timoshenko’s buckling criteria (Sigma Critical in Fig.5.7 that is found from Eq. 2.1) can be calculated. After clicking the EXECUTE button, the code runs Linear Interpolation I and Linear Interpolation II. As a result of using the code, the maximum stress at the first row of elements (MAX sigma y), total moment at the first row of elements (Moment), and the critical taper (mcr) that will result in a wrinkle for that specific element for that specific case are given as an output. 103 Figure 5.7 shows that for that element the stress calculated from Timoshenko’s buckling criteria is 247 psi. After linear interpolations, the code finds a maximum stress of 245 psi, which is very close to the Timoshenko buckling criteria. The critical taper for the case shown in Fig. 5.7 is found to be 0.00127 (in/in). The code calculated this case within six minutes. 5.5 Comparison of Results with Previous Works for the Tapered Roller Beisel [12] performed experiments and obtained data and compared his experimental results with his model. The results from the Excel VBA code will be compared with his FE model results and his experimental results. The first web he tested was a 92 gage (0.00092”) opaque polyester with a Young’s Modulus of Ex = 712000 psi, a Poisson’s Ratio of ν = 0.3, and a width of W = 6”. The nominal radius of his tapered roller was Ro = 1.49”. The compressive stress required by Eq.2.1 was approximately 265 psi for this case. Comparison of our results with his FE model and his experimental results are given below. 104 0 0.0005 0.001 0.0015 0.002 0 10 20 30 Mcr (in/in) Tension 92 ga Polyester, L=10" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.8 Wrinkles Due to Taper, 92 ga Polyester, L=10” 0 0.0005 0.001 0.0015 0.002 0 5 10 15 20 25 30 Mcr(in/in) Tension 92 ga Polyester, L=20" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.9 Wrinkles Due to Taper, 92 ga Polyester, L=20” 105 0 0.0005 0.001 0.0015 0.002 0.0025 0 5 10 15 20 25 Mcr(in/in) Tension 92 ga Polyester, L=30" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.10 Wrinkles Due to Taper, 92 ga Polyester, L=30” 0 0.0005 0.001 0.0015 0.002 0 5 10 15 20 25 Mcr(in/in) Tension 92 ga Polyester, L=40" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.11 Wrinkles Due to Taper, 92 ga Polyester, L=40” 106 Figure 5.811 shows a good agreement between our model results and his model results and his experimental results. For the 10” case some drift from the data can be observed. For longer spans (20”, 30” and 40”) our model tracks correctly with changes in span length. He also compared his model results with other test results. This was a relatively thin web, a 56 ga Polyester, with a Young’s Modulus of Ex = 658000 psi, an assumed Poisson’s ratio of ν = 0.3, and a width of W = 6”. The critical stress calculated from Eq.2.1 is 150 psi. The results from the Excel VBA code will be compared with both his model results and his experimental results below. 0 0.0005 0.001 0.0015 0.002 0 5 10 15 20 25 Mcr(in/in) Tension 56 ga Polyester, L=10" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.12 Wrinkles Due to Taper, 56 ga Polyester, L=10” 107 0.0002 0.0003 0.0008 0.0013 0.0018 0 10 20 30 Mcr(in/in) Tension 56 ga Polyester, L=20" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.13 Wrinkles Due to Taper, 56 ga Polyester, L=20” 0 0.0005 0.001 0.0015 0.002 0 5 10 15 20 25 Mcr(in/in) Tension 56 ga Polyester, L=30" Beisel's Model (COSMOS) B
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Development of Nonlinear Finite Element Codes for Stability Studies of Web Material with Strain State Dependent Properties 
Date  20111201 
Author  Yurtcu, Hasan Huseyin 
Keywords  Fem, Troughs, Vba, Web, Web Handling, Wrinkling 
Department  Mechanical Engineering 
Document Type  
Full Text Type  Open Access 
Abstract  e objective of the proposed research is to develop user friendly finite element codes that will solve nonlinear instability problems associated with strain state dependent material properties and boundary conditions of moving webs. This code will be unique and will have economic value by helping minimize web material losses. 
Note  Dissertation 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  DEVELOPMENT OF NONLINEAR FINITE ELEMENT CODES FOR STABILITY STUDIES OF WEB MATERIAL WITH STRAIN STATE DEPENDENT PROPERTIES By HASAN H.YURTCU Bachelor of Mathematical Engineering Yildiz Technical University Istanbul, Turkiye 2003 Master of Science Mathematical Engineering Yildiz Technical University Istanbul, Turkiye 2006 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of DOCTOR OF PHILOSOPHY December, 2011 ii DEVELOPMENT OF NONLINEAR FINITE ELEMENT CODES FOR STABILITY STUDIES OF WEB MATERIAL WITH STRAIN STATE DEPENDENT PROPERTIES Dissertation Approved: Dr. J.K. Good Dissertation Adviser Dr. R.P.Singh Dr .A.K.Kalkan Dr.B.Binegar Outside Committee Member Dr.Sheryl A.Tucker Dean of the Graduate College iii TABLE OF CONTENTS Page LIST OF SYMBOLS .................................................................................................................. v LIST OF FIGURES .................................................................................................................. vii CHAPTER Page 1. CHAPTER I ............................................................................................................. 1 1.1 INTRODUCTION ................................................................................................... 1 2. CHAPTER II ........................................................................................................... 7 2.1 LITERATURE REVIEW ........................................................................................ 7 2.2 Theory of Membrane Instability .............................................................................. 7 2.3 Applications of Instability Analyses to Webs ....................................................... 11 2.4 Summary ................................................................................................................ 22 2.5 Statement of Proposed Research ........................................................................... 23 3. CHAPTER III ........................................................................................................ 24 3.1 FINITE ELEMENT EQUATIONS ....................................................................... 24 3.2 Two Dimensional Four Node Quadrilateral Elements .......................................... 24 3.3 Numerical Integration by Gaussian Approach ...................................................... 33 3.4 The Global Stiffness Matrix and Matrix Banding ................................................. 37 3.5 A [D] Matrix for Plane StressIsotropic Material and Partly Wrinkled Membrane Elements ................................................................................................................ 42 4. CHAPTER IV ........................................................................................................ 50 4.1 MISALIGNED ROLLERS .................................................................................... 50 4.2 Beisel’s Method for Modeling Wrinkles Due to the Misaligned Rollers .............. 50 4.3 A New Algorithm for Predicting Wrinkles Due to the Misaligned Rollers .......... 54 4.4 Comparison of Results with Previous Works ........................................................ 66 4.5 Improvement of Execution Time ........................................................................... 68 4.6 Automating the Code ............................................................................................. 72 4.7 A New Slack Edge Criteria for a Misaligned Roller ............................................. 84 5. CHAPTER V ......................................................................................................... 91 5.1 TAPERED ROLLERS .......................................................................................... 91 5.2 Description of Tapered Rollers and Beisel’s Method for Modeling Wrinkles Due to iv CHAPTER Page the Tapered Rollers ................................................................................................ 91 5.3 A New Algorithm for Predicting Wrinkles Due to the Tapered Rollers ............... 95 5.4 Automating the Code ........................................................................................... 101 5.5 Comparison of Results with Previous Works for the Tapered Roller ................. 103 6. CHAPTER VI ...................................................................................................... 109 6.1 HOLES IN WEBS ............................................................................................... 109 6.2 Mallya’s Method for Modeling Wrinkles Due to Circular Discontinuity ........... 109 6.3 A New Algorithm for Predicting Wrinkles Due to a Circular Discontinuity ...... 113 6.4 Automating the Code ........................................................................................... 119 6.5 Comparison of Results with Previous Works for Wrinkles Due to Holes .......... 121 7. CHAPTER VII..................................................................................................... 128 7.1 NON UNIFORMITIES IN WEBS ...................................................................... 128 7.2 Modeling Non uniformities with Commercial Finite Element Code COSMOS . 129 7.3 A New Algorithm for Predicting Wrinkles Due to the Circular Non uniformities132 7.4 Comparison of COSMOS and Excel VBA Results ............................................. 138 7.5 Automating the Excel VBA Code ....................................................................... 142 8. CHAPTER VIII ................................................................................................... 145 8.1 CONVERGENCE CHECK ................................................................................. 145 8.2 The Misaligned Roller Case ................................................................................ 146 8.3 The Tapered Roller Case ..................................................................................... 154 8.4 The Central Circular Void Case .......................................................................... 161 8.5 The Central Circular Non Uniform Web Case .................................................... 166 9. CHAPTER IV ...................................................................................................... 169 9.1 SUMMARY ......................................................................................................... 169 9.2 RECOMMENDATIONS FOR THE FUTURE STUDY .................................... 170 REFERENCES ....................................................................................................................... 172 APPENDIX EXCEL VBA CODE FOR MISALIGNED ROLLERS .................................... 178 v LIST OF SYMBOLS a Web Span Length E Young’s Modulus FE Finite Elements Fig. Figure MD Machine Direction CMD Cross Machine Direction R Radius of the roller t Web Thickness Tw Web line tension X,Y,Z Coordinates σ Stress ycr σ Critical Buckling Stress md σ Machine Direction Stress cmd σ Cross Machine Direction Stress ν Poison’s Ratio ε Strain μ Coefficient of Friction I, Iz Area Moment of Inertia mcr Critical Roller Taper T Web Tension Tslack Tension Where Slack Edge Occurs vi V, Vavg Velocity, Average Velocity ω Angular Velocity σ1, σ2 Principal Normal Stresses Ffriction Friction Force FE Finite Elements x,y,z Coordinate Directions σy Normal Stress εx, εy Normal Strains ε1, ε2 Principal strains vii LIST OF FIGURES Figure Page Figure 1. 1 Troughs and Wrinkles .............................................................................................. 2 Figure 1. 2 Regime I and Regime II Wrinkles ........................................................................... 3 Figure 2. 1 Boundary Conditions for Misaligned Roller .......................................................... 11 Figure 2. 2 Webb’s Model ........................................................................................................ 15 Figure 2. 3 Webb’s Model: Boundary Conditions and Applied Loads .................................... 16 Figure 3.1 Four Node Quadrilateral Element .......................................................................... 26 Figure 3. 2 Quadrilateral Elements Inξ ,η Space (Master Element) ....................................... 27 Figure 3. 3 Two Points Gauss Quadrature using 2x2 rule. ....................................................... 36 Figure 3. 4 A Simple Model for Global Stiffness Matrix ......................................................... 37 Figure 3. 5 A Quadrilateral Element and Nodal Displacements .............................................. 38 Figure 3. 6 Mohr’s Circle for Plane Strain ............................................................................... 44 Figure 3. 7 Behaviors of Wrinkling Membrane Elements [12] ................................................ 48 Figure 4.1 Beisel’s FE Wrinkle Model for Misaligned Roller ................................................. 51 Figure 4.2 L = 6’’ Span Results ................................................................................................ 53 Figure 4.3 L = 18’’ Span Results .............................................................................................. 53 Figure 4.4 L = 30’’ Span Results .............................................................................................. 54 Figure 4.5 The System That is Modeled................................................................................... 58 Figure 4.6 Misaligned Roller FE Wrinkle Model .................................................................... 59 Figure 4.7 Flow Chart for New Algorithm ............................................................................... 61 Figure 4.8 Convergences of P and Q ........................................................................................ 63 Figure 4.9 Excel Input for the Program .................................................................................... 64 Figure 4.10 Excel Output for Stresses ...................................................................................... 65 Figure 4.11 Comparison of 6’’ Span Results ........................................................................... 67 Figure 4.12 Comparison of 18’’ ............................................................................................... 67 Figure 4.13 Comparison of 30’’ ............................................................................................... 68 viii Figure Page Figure 4.14 Flow Chart of Computer Code for Improved Execution Speed ............................ 71 Figure 4.15 Linear Interpolation I ............................................................................................ 74 Figure 4.16 Flow Chart for Linear Interpolation I.................................................................... 76 Figure 4.17 Linear Interpolation II ........................................................................................... 78 Figure 4.18 Flow Chart for Linear Interpolation II .................................................................. 79 Figure 4.19 Excel Input for the Automated Program ............................................................... 81 Figure 4.20 Comparison of 6’’ Span Case with the Modified Code ........................................ 82 Figure 4.21 Comparison of 18’’ Span Case with the Modified Code ...................................... 82 Figure 4.22 Comparison of 30’’ Span Case with the Modified Code ...................................... 83 Figure 4.23 Slack Edge ............................................................................................................. 84 Figure 4.24 Slack Edge Criteria ............................................................................................... 85 Figure 4.25 Beams with Shear .................................................................................................. 86 Figure 4.26 Moment in a Beam ................................................................................................ 87 Figure 5.1 Tapered Roller Profile ............................................................................................. 92 Figure 5.2 Beisel Tapered Roller Model [12] .......................................................................... 93 Figure 5.3 The Figure for Tapered Roller ................................................................................ 96 Figure 5.4 Tapered Roller FE Wrinkle Model ......................................................................... 96 Figure 5.5 The Flow Chart for Tapered Roller ......................................................................... 99 Figure 5.6 Calculating Moment and Critical Taper ................................................................ 100 Figure 5.7 Excel Input of Tapered Roller ............................................................................... 102 Figure 5.8 Wrinkles Due to Taper, 92 ga Polyester, L=10” ................................................... 104 Figure 5.9 Wrinkles Due to Taper, 92 ga Polyester, L=20” ................................................... 104 Figure 5.10 Wrinkles Due to Taper, 92 ga Polyester, L=30” ................................................. 105 Figure 5.11 Wrinkles Due to Taper, 92 ga Polyester, L=40” ................................................. 105 Figure 5.12 Wrinkles Due to Taper, 56 ga Polyester, L=10” ................................................. 106 Figure 5.13 Wrinkles Due to Taper, 56 ga Polyester, L=20” ................................................. 107 Figure 5.14 Wrinkles Due to Taper, 56 ga Polyester, L=30” ................................................. 107 Figure 5.15 Wrinkles Due to Taper, 56 ga Polyester, L=40” ................................................. 108 Figure 6.1 Troughs and wrinkles due to a circular hole in the web [41] ................................ 110 Figure 6.2 Distance Between the Hole and the Tangent Line ................................................ 111 Figure 6.3 Mallya’s Model ..................................................................................................... 112 ix Figure Page Figure 6.4 The Figure for the Hole Model Traveling Between Rollers ................................. 113 Figure 6.5 The Web Hole FE Wrinkle Model ........................................................................ 114 Figure 6.6 Meshing the Hole Region...................................................................................... 115 Figure 6.7 Excel Output of Hole Region ................................................................................ 116 Figure 6.8 The Flowchart for Circular Discontinuity ............................................................. 118 Figure 6.9 The Section of Web over the roller ....................................................................... 120 Figure 6.10 Comparison of Results for Hole .......................................................................... 121 Figure 6.11 Excel output of square ......................................................................................... 122 Figure 6.12 Excel output of equilateral quadrangle ............................................................... 122 Figure 6.13 Comparison of Circle, Equilateral Quadrangle and Square Shaped Holes ......... 123 Figure 6.14 Flow Chart for Automating Hole Code ............................................................... 125 Figure 6.15 Excel Interface of Circular Void Excel VBA Code ............................................ 127 Figure 7. 1 Non uniform Web Models (COSMOS) ............................................................... 130 Figure 7. 2 Critical CMD Stresses Developed When the Non uniform Area is 3” Away From the Roller............................................................................................................... 131 Figure 7. 3 A Web with a Non uniform Circular Shaped Region Travelling Between Rollers132 Figure 7. 4 The Model for the Non uniform Hole Shaped Material ....................................... 133 Figure 7.5 Meshing the Non uniform Region......................................................................... 134 Figure 7. 6 Troughs and Wrinkles .......................................................................................... 135 Figure 7. 7 The Flowchart for Non uniform Webs ................................................................. 137 Figure 7. 8 Wrinkles Due to Non uniformities, 2r = 3”, t’ = 0.0008” .................................... 139 Figure 7. 9 Wrinkles Due to Non uniformities, 2r = 2”, t’ = 0.0008” .................................... 139 Figure 7. 10 Wrinkles Due to Non uniformities, 2r = 3”, t’ = 0.0001” .................................. 140 Figure 7. 11 Wrinkles Due to Non uniformities, 2r = 2”, t’ = 0.0001” .................................. 140 Figure 7. 12 Wrinkles Due to Non uniformities, L = 3”, t’ = 0.0005” ................................... 141 Figure 7. 13 Excel Interface of Non uniform Excel VBA Code ............................................ 143 Figure 8. 1 Meshing the Model .............................................................................................. 147 Figure 8. 2 The Flowchart for the Convergence Check ......................................................... 149 Figure 8. 3 Comparison of 6” Span Case after the Convergence Check ................................ 150 Figure 8. 4 Comparison of 18” Span Case after the Convergence Check .............................. 151 Figure 8. 5 Comparison of 30” Span Case after the Convergence Check .............................. 151 x Figure Page Figure 8. 6 Converging Result for a Specific Case for Misaligned Roller Case .................... 152 Figure 8. 7 Converging Result for a Long Narrow Web ........................................................ 153 Figure 8. 8 Converging Result for a Wide Web ..................................................................... 153 Figure 8. 9 Converging Result for 0.002” Thick 50” Wide Web ........................................... 154 Figure 8. 10 The Flowchart for the Tapered Roller ................................................................ 156 Figure 8. 11 Converging Results for 20” Case ....................................................................... 157 Figure 8. 11 Converging Results for 30” Case ....................................................................... 158 Figure 8. 11 Converging Results for 40” Case ....................................................................... 158 Figure 8. 14 Converging Result for a Specific Case for Tapered Roller Problem ................. 159 Figure 8.15 Converging Result for a 60” Long Web ............................................................. 160 Figure 8.16 Converging Result for a 80” Long Web ............................................................. 160 Figure 8. 17 Convergence Check for the Void Case .............................................................. 163 Figure 8. 18 Result of Modified Code for the Case at Figure 6.15 ........................................ 164 Figure 8.19 Convergence Result for a Wide Web for Circular Void Code............................ 165 Figure 8.20 Convergence Result for a Long Web for Circular Void Code ............................ 165 Figure 8. 21 Result of Modified Code for the Case at Figure 7.15 ........................................ 167 Figure 8. 22 Result of Modified Code for a Wide, Long Web ............................................... 168 1. CHAPTER I 1.1 INTRODUCTION A web is any flexible thin material. Webs are made in continuous production processes. The webs are often stored in wound roll form, since this is the only convenient means of storing long lengths of flexible material. A production roll of polyester film that is 48 gages (0.00048”) in thickness might be 4 feet in diameter. If the web is wound on a 6” diameter core there can be nearly 60 miles of web length in the roll. That polyester web could be nearly 300 inches wide when made, but costumer rolls to be shipped could be only 6 inches to 60 inches in width depending on the products to be made. Thus a particular web might be 60 miles long, 60 inches wide, and 0.00048 inches in thickness. Plastic films, papers, foils, and thin metal sheets are examples of webs. Web handling can be defined as all processes employed during the transportation of webs. Cutting, coating, slitting, printing, laminating and drying are some of the processes that add value to the web. During these processes, web materials must travel around several rollers in the process equipment. 2 During the transportation of webs trough process machinery, compressive stresses are induced in webs and these compressive stresses may cause instabilities. If these instabilities occur in free spans, which are the sections of web between two rollers, they are called troughs and if they occur at rollers, they are called wrinkles (Figure 1). Wrinkles and troughs result in loss of value and quality of the webs. For instance, if wrinkles may cause permanent damages to the webs the result is wasted web and highly priced downtime of web lines. The processes such as laminating, printing on paper, metalizing of films or coating require web to be planar. Troughs may affect these processes and like wrinkles the result is waste of material and time. As a result, wrinkles and troughs are two important engineering problems in the web handling industry. Figure 1. 1 Troughs and Wrinkles There are three types of wrinkles. Machine direction wrinkles, cross machine direction wrinkles and shear wrinkles. The machine direction [MD] is the direction of the web travel through the web process machine and the cross machine direction [CMD] is the direction that is perpendicular to the machine direction (i.e. across the web width). MD wrinkles occur 3 because of the compressive forces in the CMD. CMD wrinkles are produced during winding or unwinding processes due to interlayer slippage that may be the result of air entrainment [1]. Shear forces due to the imperfect rollers, misaligned rollers and non uniform webs may result in shear troughs and wrinkles. Shear wrinkles can result from misaligned or tapered rollers [1, 2]. Shear wrinkles can be classified as regime one wrinkles and regime two wrinkles. Roller Misalignment or Taper Web Tension Regime I Regime II Planar, Troughed, or Slack Edge Spans Wrinkles Figure 1. 2 Regime I and Regime II Wrinkles Regime one wrinkles are diagonal wrinkles which occur because of the presence of a lateral shear force in the web. The shear force might be due to a misaligned or tapered roller. The 4 shear forces in a web which are necessary to cause troughs or wrinkles are very small. The misalignments and tapers which cause the shear forces are unintentional. For regime one troughs and wrinkles to occur there must be sufficient traction between the web and the rollers to react the shear forces. The shear forces occur because when traction is adequate the web will attempt to gain normal entry to the downstream roller. The concept or law of normal entry of a web to a roller is attributed to Lorig [3]. Regime two wrinkles are dependent on traction and velocity between the web and the roller. In this case, if traction is insufficient, CMD compressive stresses cannot develop in the web. Frictional restraint forces between the web and the roller are required to sustain wrinkles. Increases in web tension result in increased normal forces between the web and roller and act to decrease the air films which develop between web and roller. Air films can develop between webs and rollers due to the hydrodynamic entrainment of air. Increased web tension and normal force act to decrease the thickness of the air films. The air films result in decrease in traction between the web and rollers. Thus an increase in web tension will increase the potential frictional forces that would be needed to sustain a wrinkle on a roller. The normal entry rule may be violated as a result [3]. Good et al [4] showed that regime two wrinkles could be predicted. In this proposal it will be assumed that sufficient traction exist between webs and rollers such that only regime one wrinkles need to be considered. Wrinkles are defined as buckling of webs around rollers. Webs have very small resistance to compression in free spans. Webs may withstand more compressive stress on the rollers than free span, because the web has greater stability in the form of a cylindrical shell than it has a flat plate. The critical CMD compressive stress needed to wrinkle the web can be predicted using classical cylindrical shell buckling expressions. 5 Previous research by Webb[9], Beisel[12] and Mallya[15] have shown that web wrinkles on rollers can be predicted by analyzing the compressive CMD stresses that form in the web as it contacts rollers. These compressive stresses arise due to troughs that have already formed in the free span prior to the rollers. The troughs themselves may be the result of many disturbances that exist in web lines. Misaligned rollers, tapered rollers and web non uniformities (such as holes) are examples of such disturbances. Thus, it could be stated that web wrinkles on rollers are the result of two instabilities. A disturbance is first required to induce a trough instability in the free span. The first instability has now occurred. The web that was planar in the free span has now troughed. After the troughs appear CMD compressive stresses will arise in the web entering the roller downstream of that span. When the troughs first appear, the CMD stress in the web entering the roller may be small. Thus, whatever disturbance produced the trough initially may have to become yet larger before sufficient CMD compressive stress can be generated in the web on the roller to produce an MD wrinkle. At this point, the second instability has occurred. The web on the roller that earlier had the shape of a sector of a cylindrical shell has now buckled. Thus, it has been proven that the prediction of wrinkles upon a roller involves: (1) An instability (troughs) occurring in the span upstream of the roller. (2) As whatever disturbance increases that produced the troughs in step (1), a post buckling analysis must be undertaken. It will bill seen that as the disturbance increases that a CMD compressive stress will arise in the web on the downstream roller. (3) As the disturbance yet increase further the CMD compressive stress in the web on the 6 downstream roller is also increasing. As this CMD stress increases it will finally surpass the cylindrical shell buckling stress. At this point, with sufficient friction between the web and the roller, a wrinkle will form in the web on the roller. The purpose of this study is to develop efficient computational tools that can accomplish the analyses required in step (1), (2), and (3). 7 2. CHAPTER II 2.1 LITERATURE REVIEW The literature has been reviewed and the findings will be broken into two sections. First the basic theories of membrane instability will be reviewed. Second, those studies which examine web instability will be reviewed. At the close of this chapter a final summary of the findings will be included and a statement of proposed research will be presented. 2.2 Theory of Membrane Instability Wagner [18] prepared a treatise on sheet metal girders with very thin webs. Probably this study is the earliest investigation of the mechanics of wrinkling membranes. He worked to develop the structural method of sheet metal girders. His methods were based on the 8 assumption of the low stiffness in bending of the metal web. He worked to explain the behavior of the thin metal webs in beams carrying a shear load well in excess of the initial buckling value. He proposed tension field theory. Tension field theory helps to analyze flexible structures that can support tension, but have no ability to resist compression. This theory was further developed by Reissner [23], SteinHedgepeth [19] and Mikolas [20, 21]. MillerHedgepeth [5] and Miller et al. [6, 22] adopted this theory to the finite element method. Tension field theory can be applied to the web lines due to the fact that web lines can support tension but cannot carry compression and also web materials are flexible structures. Stein and Hedgepeth [19] suggested a particularly useful approach concerning partly wrinkling membranes. This work is a seminal work in this field. They derived a theory to predict the stresses and deformations of stretched membrane structural components for loads under which part of the membrane wrinkles. Their theory was based on the basic assumption that a membrane has no bending stiffness and because of this can carry no compressive stress. They applied their theory to inplane bending of a stretched rectangular membrane, a pressurized cylinder, and to the rotation of a hub in a stretched circular membrane. They presented stresses and deformations in equation form for the wrinkled and unwrinkled regions. The membrane they considered is elastic, isotropic, has no bending stiffness, and cannot carry compressive stress. In their work they considered average strains and displacements of the wrinkling material rather than detailed deformations of each wrinkle. In terms of the wrinkling equations given, their theory was limited in the sense that average strains must be small compared with unity. They started to investigate the wrinkling region 9 by looking the principal stresses. They used a criterion that if both principal stresses are zero, the membrane is unloaded and thus will not wrinkle. Their criteria for a wrinkled membrane is one principal stress must be zero and other principal stress is non zero and tensile. The nonzero principal stress may be assumed to act along the crest of the wrinkle. The approach that was developed by Stein and Hedgepeth [19] was further developed and applied by Mikolas [20, 21]. He presented experiments and analysis for the wrinkling behavior of stretched membranes under the influence of a torque loading through an attached hub. He found that theory and experimental results were in a very good agreement. The work done by Stein and Hedgepeth [19] and Mikolas [20, 21] were closed form solutions and did not involve the finite element method. The principal stress criterion that they used is employed in our current approach. Miller and Hedgepeth [5] developed a new algorithm for finite element analysis based on the same assumptions and field equations after finding some critical disadvantages connected with the Stein and Hedgepeth [19] approach. This work may be the most important study in this field. In their algorithm the element stiffness is dependent on the current principal strains. Wrinkling membrane elements can have either taut behavior, wrinkled behavior or slack behavior. In taut behavior both principal stresses are larger than zero, in wrinkled behavior one of the principal stresses is greater than zero the other is equal to zero and in slack behavior both principal stresses are equal to zero. In their algorithm, in the first load step all elements are assumed to taut behavior. In the consecutive steps, element behaviors 10 are calculated with respect to strain states of previous steps. In other words, the decision on the stress state is made using the criteria based upon principal strains. In the algorithm they apply load step by step and they continue to solve for a particular load step until the convergence is achieved for that load step. In the proposed research, the algorithm that is employed is similar to Miller and Hedgepeth’s algorithm. Miller et al. [6, 22] investigated the algorithm further. They presented the efficiency and accuracy of the algorithm by applying it to the problems Stein and Hedgepeth [19] studied. They described the algorithm more detailed in these studies. They also described how Miller and Hedgepeth [5] derived the Dtaut, Dslack and Dwrinkled constitutive matrices in detail. 11 2.3 Applications of Instability Analyses to Webs Shelton [17] studied the steering effects of downstream misaligned roller. He modeled web span as a beam. His work helped us to justify boundary conditions which should be associated with a straight uniform web approaching a misaligned roller when there is sufficient friction to enforce the kinematic boundary conditions. Figure 2.1 Boundary Conditions for Misaligned Roller He determined that the concept of normal entry that had been used in the drive belt industry can also be used in web guiding systems. He stated that Vi and θi can be arbitrary set equal to zero without effecting the relative lateral deformation within the span. θj is the misalignment of the downstream roller at j, and due to the law of normal entry the slope of the web at j will be θj. Shelton proved that the final boundary condition was the downstream moment being zero. Thus Mj(L) = Vj’’(L) = 0. In our problem these boundary conditions will be used while modeling the instability of a web due to a downstream 12 misaligned roller. Gehlbach, et al. [32] proposed buckling criteria for troughed webs in a free span. They showed experimental verification for downstream misaligned roller case. In their work, isotropic web properties were considered. Gopal and Kedl [33] were the fist people who used finite element analysis and a commercial FE code to study trough formation in the web span between rollers. They modeled a web span by using triangular plate elements and they used ABAQUS FE commercial code. They were successful to predict deformations of a web due to the misaligned roller. Benson et al. [8] developed a finite element model for wrinkling analysis and called this code FEWA. They used this code to make calculations in their paper. Their aim was to better understand the conditions which cause wrinkle formation. They worked to predict locations where troughs would form and predict magnitudes of compressive stresses. They compared some of their results with the results from the nonlinear version of commercial code ABAQUS. In their code instead of using Wagner’s tension field theory they used the tension field theory that Wu [24, 25, and 26] introduced first. In this method, it is assumed that out of plane deflection relieves compressive stress across a wrinkle and that there is an associated strain with this deflection. Roisum [28, 29, 30, and 31] described wrinkling phenomena in detail. He explained 13 wrinkling, air entrainment, tension control, roller design, problems associated with profile variations, why and how wrinkles form, the types of wrinkles, troubleshooting techniques. He explained the importance of the problem due to the material cost and waste for the producers. While doing this he explained the subject in a simple way rather than an academic way. Shelton [14] worked on buckling of webs. He modeled a web as a buckled plate or a shell. He studied buckled wavelengths of webs in free spans and on rollers. He predicted the wavelength of the buckled form and he compared these results with experimental data. Good et al. [4] worked on velocity independent and velocity dependent wrinkles. Velocity dependence occurs when velocities are high enough and web tension is low enough to allow sufficient air entrainment between webs and rollers that Shelton’s boundary conditions [17] are no longer valid. They found that velocity dependent wrinkles can be avoided by using enough traction and suitable web line conditions. Hashimoto [34] worked on the studies done by Beisel [2, 10], which is a theoretical model for predicting the web wrinkling due to the misaligned roller. The theoretical model was established upon the experimental results. The experiments which he did were for non uniform webs with different Young’s modulus in MD and CMD directions. He compared these experimental results with the model and he verified the accuracy of the model. 14 Jones and McCann [35] studied wrinkling of webs on rollers and drums. They used a new variational analysis based on the method of RayleighRitz. They modeled the wrinkling as a continuous sine wave in a web on a roller or a drum. They reported that a shell buckling on a rigid support in an outward mode is similar to wrinkling on rollers. Beisel [2, 10] studied wrinkling phenomena due to the misaligned roller. In his study, he considered orthotropic material properties. He showed that while wrinkle formation would require post buckling analysis, prediction of trough formation can be expressed by closed form solution. Papandreadis [7] employed finite element methods to predict troughs in the webs. He studied effects of several parameters on the amount of lateral contraction of the web. He examined the effects of web material properties (Poisson’s ratio, modulus of elasticity) and web geometry (various lengthtowidth ratios), web thickness, loading conditions (tensile loading, combination of tension and shear forces) on the wrinkling phenomenon. He used a finite element code, named NASTRAN (Nasa Structural Analysis), to analyze the buckling of panels and the resulting shapes (like wavelength of the corrugations, number of waves). Webb [9] was the first person who tried to couple the behavior of the web in the free span to the web on the roller. Probably, the most important finding in his work was that CMD compressive stresses were forming in the web on the roller due to the troughs that had already formed in the free spans. He used quadrilateral elements within the commercial finite element code COSMOS to predict wrinkles due to the downstream misaligned roller. 15 He observed that there is a linear relationship between wrinkle θ and trough θ . He performed experiments to measure the deflection as a result of a misaligned roller. He found that the misalignment required to form a wrinkle was roughly twice that to form a trough over a wide range of parameters. He tested a wide range of web material properties, web thickness, and span length to width ratio (often called span ratio). He tried to find a relationship between the width reduction of the web and web buckling. He modeled the web crossing the roller using regular elastic elements. He used wrinkling membrane elements for the web in free span. In Fig.2.2 Webb’s approach to the problem is presented. Figure 2.2 Webb’s Model In Fig.2.3 the boundary conditions and the applied loads in Webb’s model is shown. 16 Force Prescribed Displacement Coupled DOFs Fixed DOF X Y Figure 2.3 Webb’s Model: Boundary Conditions and Applied Loads Timoshenko [11] showed the axial buckling stress for a sector of a cylinder is the same per unit circumferential length as that of a complete cylinder. He found that ycr σ is 3 (1 2 ) ycr h E R σ ν = − ⋅ ⋅ − (2.1) The web wrapped around the roller can be assumed to be sector of a cylindrical shell. Shelton [14] studied with the mechanics of buckling and he hypothesized that the web wrapped around the roller is shows similarity to an internally pressured thin cylinder vessel. He found that the tension in the web performs like the hoop stress in a pressurized cylinder and that a 17 pressure results between the web and roller. He discussed that Eq.2.1 may indeed be appropriate to use in modeling a web buckling on a roller. Good and Beisel [10] developed an instability criterion for orthotropic web in shell form in order to predict wrinkling of webs: 3(1 ) x y ycr x y t E E R v y σ = − − (2.2) where R is the radius of the cylinder, t is the thickness of the web, Ex and Ey are the elastic modulus in x and y directions respectively, and v is Poisson’s ratio. Beisel performed many tests of aluminum soda cans which were near perfect in geometry. He found that the buckling stress approached that given by Timoshenko’s expression and concluded that the earlier disagreement being due to shell imperfections must have been correct. Thus expressions 2.1 and 2.2 appear to be applicable to sectors of web transiting rollers. If the compressive stresses in the web on the roller reach the ycr σ , the web will buckle. Webb [9] increased the shear forces used to simulate the misalignment of the roller until the compressive stress induced in the web on the roller reaches the ycr σ value in (Eq.2.2). When the critical compressive CMD stress was reached, he determined the rotation of the downstream roller from the displacements output by the finite element code. The experimental critical misalignment of the roller and the rotation of the web at the entry of the 18 roller computed using finite element model did not match well. He proposed that the experimental critical misalignment of the roller cr θ is the sum of the critical angles predicted by the model for troughing, cr ,trough θ and for wrinkling, cr ,wrinkle θ in verifying his model. This was found later to be incorrect by Beisel [12]. Webb’s work was carried forward by Beisel [12]. Beisel studied troughing and wrinkling due to roller misalignment, roller taper, and roller crown. There is some similarity between the effects of roller misalignment and taper on troughs and wrinkles. In both cases the web is steered laterally in the machine and a shear stress results in the web. The misaligned roller induces the misalignment angle and a lateral deformation at the misaligned roller. The tapered roller induces a lateral deformation and a bending moment at the tapered roller. Beisel’s method of modeling the problem with COSMOS was different from the Webb’s method of modeling the problem. Beisel achieved good agreement with his models and the experimental results. His model and his results will be introduced in Chapter IV and Chapter V. Swift [39] examined steering of drive belts. He worked with the concept of a couple developed in a web approaching a tapered or crowned roller and the resulting steering of the web. He suggests the minimum amount of taper or crown which should be employed to control the web with minimum interference of stresses in the web. He gave experimental results to support his suggestions for corrective measures. 19 Shelton [44] discussed misaligned and tapered rollers, and he mentioned multiple web span interactions and moment transfer. In his thesis, Beisel [12] compared his model with Shelton’s model for the critical taper required to cause troughs to form in a web span approaching a tapered roller. Good and Beisel [45] worked on the formation of troughs due to tapered rollers. They attempted to determine whether the procedure that was employed for misaligned rollers would be applicable to the case of downstream tapered rollers. This work was a part of Beisel’s thesis [12]. Brown [46] presented a new method for modeling the elastic behavior of webs conveyed over rollers. He worked on lateral displacement of a misaligned roller and lateral displacement of a tapered roller. He suggested two modifications of the web boundary conditions. One was a generalization of the normal entry rule and the other was the addition of what he named the normal strain rule. With a numerical partial differential equation solver, he solved two dimensional plane stress equations and compared the results with earlier models. Shimizu et al. [47] and Shimizu [48] worked on plates which have holes and are subjected to tensile load. They used the finite element method in their work. They investigated the effects of aspect ratios (height/width) and shapes of holes to the k which is the buckling coefficient 20 of the plate described by Timoshenko and Gere [11]. They determined that with that curvatures on corners of holes have little effect in improving the tension buckling strength and that the buckling coefficient increases corresponding to the increasing aspect ratio. ElSawy and Martini [49] studied the effects of plate aspect ratio (height/width), hole size, hole location and loading ratio on the buckling coefficient k of rectangular perforated plates subjected to uniform end compression in the x direction and compression or tension in the y direction. Mallya [15] was the first person to examine the effects of holes in process machinery in web handling industry. He applied Beisel’s method of modeling to web wrinkles due to circular and elliptical discontinuities in the web. He studied the behavior of webs with voids traveling over a roller. He compared experimental results and finite element model results that he modeled using commercial FE code COSMOS. He studied elliptical voids and circular holes in terms of generating wrinkles. His FE model was similar in form to Beisel’s FE model with respect to boundary conditions and using five panels. His model and results will be introduced in Chapter VI. Kara [16] also used similar modeling method to predict the occurrence of wrinkles due to length variation across the web width. He attempted to find the critical conditions that would induce wrinkles. He heated the center of the web during his experiments to achieve length variation. He also used COSMOS FE commercial code to model his case, and used experimental findings to confirm his FE model. 21 Relevant books include, ‘Theory of Elastic Stability’ written by Timoshenko and Gere [11], ‘Introduction to Finite Elements in Engineering’ written by Chandrupatla and Belegundu [36], ‘Finite Elements and Approximation’ written by Zeinkiewicz and Morgan [37] and Visual basic Excel for Dummies [27]. Also ‘Finite Element Analysis Class Notes’ from Good [38] is used for finite element part of the study. 22 2.4 Summary The tension field theory can be applied to webs in web lines because the web can support tension but cannot carry compression. Wagner [18] proposed tension field theory. His theory was further developed by Raissner [23], Stein and Hedgepeth [19] and Mikolas [21, 22]. Miller and Hedgepeth [5] developed a new algorithm for finite element analysis based on this field equation. This algorithm is the most common algorithm used to examine wrinkles. Lorig [3] developed the normal entry rule for a web approaching a roller. Shelton [17] found that the moment in a web is zero when it approaches a misaligned roller. In the web handling area, Gopal and Kedl [33] were the first to employ a commercial finite element code to study trough formation in the web span between the rollers. Webb [9] attempted to model wrinkle formation due to a misaligned roller. He used partly wrinkling membrane elements while modeling the web span. He used the commercial finite element code called COSMOS while modeling his work, and shell buckling criteria to determine whether wrinkling occurred. Beisel [12] made the most recent attempts to model wrinkle formation due to the misaligned roller and the tapered roller. He modeled the web between the rollers by wrinkling membrane elements and linked it to the classical shell buckling criteria as Webb [9] did. He studied webs approaching misaligned rollers and tapered rollers. He used the commercial finite element code COSMOS to model these cases. He compared his model results with his experimental results .These results showed good agreement. Mallya [15] was the first person who investigated the effects of voids on the stability of 23 webs. He performed experiments and modeled circular and elliptical discontinuities in webs in web lines. He also used the commercial finite element code COSMOS in his models and his model results and test results matched well. 2.5 Statement of Proposed Research In the literature no citations were found for the wrinkling instability of moving webs that did not involve the use of commercial finite element codes for solution. Use of commercial finite element codes by novice users to solve nonlinear problems associated with web instability is difficult. The objective of the proposed research will be to develop user friendly finite element codes that will solve nonlinear instability problems associated with strain state dependent material properties and boundary conditions of moving webs. This code will be unique and will have economic value by helping minimize web material losses as described in the introduction. 24 3. CHAPTER III 3.1 FINITE ELEMENT EQUATIONS In this chapter, the finite element equations will be described in detail. Displacement, strain and stress equations, the element stiffness matrix, meshing, banded matrix and their relations will be studied. A discussion of the solution method for cases in which the elasticity matrix [D] is not constant will be given. 3.2 Two Dimensional Four Node Quadrilateral Elements In this study, two dimensional four node quadrilateral elements are used. In this section the equations and properties of four node quadrilateral elements will be given briefly. 25 The stiffness matrix for quadrilateral elements can be found from the strain energy term in a continuum. 1 2 T V U = ∫ σ ε dV (3.1) In finite elements we consider the total strain energy to be sum of the strain energies from each element. Eq.3.2 can be obtained by replacing dV by t dA in Eq.3.1, where t is the uniform thickness of an element. 1 2 T e e e U =Σt ∫ σ ε dA (3.2) The small strain displacement relations for two dimensional problems can be written as: x y xy u x v y u v y x ε ε ε γ ∂ ∂ ∂ = = ∂ ∂ ∂ + ∂ ∂ (3.3) Where u and v are the deformations in the x and y directions respectively. In twodimensional fields, the displacement components at each point in the domain of the finite element can be represented as functions of two coordinate directions. u = [u(x, y), v(x, y)]T (3.4) 26 For the general quadrilateral element shown in Fig.3.1, the nodal displacement vector is q = [q1,q2 ,q3 ,...,q8 ]T X Y u P(x,y) 4 2 1 3 q1 q2 q3 q4 q6 q8 q5 q7 v Figure 3.1 Four Node Quadrilateral Element The finite element method uses concept of shape functions to develop interpolations systematically. According to the concept, the shape functions must be developed for the master element. The master element is defined in natural coordinates ( ξ ,η ) and has a square shape in the natural coordinate system (Fig.3.2). The Lagrange shape functions are 1 2 3 4 N , N , N and N . Shape functions take the value of unity at the node where they are defined and the value zero at the 27 other nodes. For example 1 N takes the value one of unity at node one and takes the value zero at the other nodes. ξ η Figure 3. 2 Quadrilateral Elements Inξ ,η Space (Master Element) At the edges ξ = +1 andη = +1 1 N is equal to zero. So, 1 N must be a function like, 1 N = c (1−ξ ) (1−η ) (3.5) c is a constant that can be determined easily. Since 1 N is equal to one at node one, where ξ = −1 andη = −1.If we put these values at Eq.3.5, 1 = c (1− (−1)) (1− (−1)) (3.6) yields c =1/ 4 .Finally 1 N can be written as, 1 1 (1 ) (1 ) 4 N = −ξ −η (3.7) 28 By using the same procedure N1, N2 , N3 and N4 can be written as, 1 2 3 4 1 (1 ) (1 ) 4 1 (1 ) (1 ) 4 1 (1 ) (1 ) 4 1 (1 ) (1 ) 4 N N N N ξ η ξ η ξ η ξ η = − − = + − = + + = − + (3.8) Now, the shape functions can be used to interpolate the displacement at any point within the domain of the element using the equations: 1 1 2 3 3 5 4 7 1 2 2 4 3 6 4 8 u N q N q N q N q v N q N q N q N q = + + + = + + + (3.9) if N is, 1 2 3 4 1 2 3 4 0 0 0 0 0 0 0 0 N N N N N N N N (3.10) then the displacement can be written as: [ ]{ } 1 2 3 1 2 3 4 4 1 2 3 4 5 6 7 8 0 0 0 0 [u] 0 0 0 0 q q q u N N N N q N q v N N N N q q q q = = = (3.11) with the help of shape functions a point in the element can be described as, 29 1 1 2 2 3 3 4 4 1 1 2 2 3 3 4 4 x N x N x N x N x y N y N y N y N y = + + + = + + + (3.12) The shape functions can also be used to generate a map between the cartesian coordinates (x,y) and the natural coordinates (ξ ,η ). Since the same shape functions have been used to interpolate deformation within an element and to generate the coordinate map equations this is called an isoparametric formulation. The strain relationships (3.3) require derivatives with respect to cartesian coordinates. We currently have the deformations u and v defined as shape functions, which are functions of the natural coordinates ξ andη , multiplied by nodal deformations that are constants. So, to determine the strains in cartesian coordinates we must first relate the derivatives of deformations in natural coordinates (ξ ,η ) to derivatives in cartesian coordinates(x,y). If a displacement function in x , y coordinates is u = u(x , y) then this function can be considered to be an implicit function of ξ andη asu = u[x(ξ ,η ) , y(ξ ,η ) ] . Differentiation due to the chain rule, u u x u y x y u u x u y x y ξ ξ ξ η η η ∂ ∂ ∂ ∂ ∂ = + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = + ∂ ∂ ∂ ∂ ∂ (3.13) If we define the Jacobian matrix as, J x y x y ξ ξ η η ∂ ∂ ∂ ∂ = ∂ ∂ ∂ ∂ (3.14) We can rewrite Eq.3.13 as, 30 =J u x y u u x x u x y u u y y ξ ξ ξ η η η ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ (3.15) It seems that the Jacobian can transform derivatives in cartesian coordinates to derivatives with respect to natural coordinates. By using Eq.3.8, Eq.3.12 can be written as, 1 2 3 4 1 2 3 4 1 1 1 1 (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) 4 4 4 4 1 1 1 1 (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) 4 4 4 4 x x x x x y y y y y ξ η ξ η ξ η ξ η ξ η ξ η ξ η ξ η = − − + + − + + + + − + = − − + + − + + + + − + (3.16) With the help of Eq.3.16 Jacobian term Eq.3.14 can be written as, 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) J 4 (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) x x x x y y y y x x x x y y y y η η η η η η η η ξ ξ ξ ξ ξ ξ ξ ξ − − + − + + − + − − + − + + − + = − − − + + + − − − − − + + + + − (3.17) this equation can be written as, 11 12 21 22 J J J J J = (3.18) Now the Jacobian can be inverted and rewrite Eq.3.15 to produce derivatives that are related to strains: J 1 u u x u u y ξ η − ∂ ∂ ∂ ∂ = ∂ ∂ ∂ ∂ (3.19) 31 by using Eq.3.17, Eq.3.19 can be written as 22 12 21 11 1 det J u u x J J u J J u y ξ η ∂ ∂ ∂ − ∂ = ∂ − ∂ ∂ ∂ (3.20) by following the same procedure v displacements can be obtained as, 22 12 21 11 1 det J v v x J J v J J v y ξ η ∂ ∂ ∂ − ∂ = ∂ − ∂ ∂ ∂ (3.21) The equation ∂x ∂y = det J ∂ξ ∂η has a proof in reference [36]. By using Eq.3.3, Eq.3.20, Eq.3.21 and defining A matrice as, 22 12 21 11 22 11 22 12 0 0 1 0 0 det J J J A J J J J J J − = − − − (3.22) The strain displacement relations can be written as 22 12 21 11 22 11 22 12 0 0 1 0 0 det J x y xy u u u x u u J J v J J A y v v J J J J u v y x v v ξ ξ ε η η ε ε γ ξ ξ η η ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ − ∂ ∂ ∂ = = = − = ∂ ∂ ∂ − − ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ ∂ ∂ (3.23) with the help of 3.17 G matrice can be defined as, 32 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 1 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 4 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) G η η η η ξ ξ ξ ξ η η η η ξ ξ ξ ξ − − − + − + − − − + + − = − − − + − + − − − + + − (3.24) and using the displacement vector Eq.3.9 the derivatives of u and v can be written in the natural coordinates as a function of the nodal deformations, (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 1 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 4 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) u u q G q v v ξ η η η η η ξ ξ ξ ξ η η η η ξ ξ ξ ξ ξ η ∂ ∂ ∂ − − − + − + ∂ − − − + + − = = ∂ − − − + − + ∂ − − − + + − ∂ ∂ (3.25) If B is defined as B=AG, by using Eq.3.23 and Eq.3.25 strain term can be written as, 1 2 3 22 12 4 21 11 5 22 11 22 12 6 7 8 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 0 0 1 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 0 0 4det J 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) q q q J J q J J q J J J J q q q η η η η ξ ξ ξ ξ ε η η η η ξ ξ ξ ξ − − − + − + − − − − + + − = − − − − + − + − − − − − + + − B q = (3.26) two dimensional constitutive relations will be used to relate stress to strain σ =Dε and now, σ =D Bq (3.27) 33 By using Eq.3.26 and Eq.3.27 the strain energy term 1 2 T e e e U =Σt ∫ σ ε dA , becomes, 1 1 1 1 1 [ det ] 2 1 2 T T e e T e e U q t B DB J d d q q k q ξ η − − = = Σ ∫ ∫ Σ (3.28) where k is 8x8 element stiffness matrix: 1 1 1 1 e T det e k t B DB J dξ dη − − = ∫ ∫ (3.29) These integrals can be evaluated by using numerical integration methods. The Gaussian approach will be considered for this purpose. 3.3 Numerical Integration by Gaussian Approach By integrating in natural coordinates the bounds of integration are much simplified. In cartesian coordinates the ybounds will be functions of x and the xbounds will be functions of y. In natural coordinates our bounds are from 1 to1 for both ξ andη . Series can be used to take the integrals (Eq.3.29).By using the Gaussian quadrature approach; integration can be evaluated using weights and points. These points are also called Gauss points. 1 1 1 2 2 1 ( ) ( ) ( ) ( ) n n I f ξ dξ ω f ξ ω f ξ ω f ξ − = ∫ ≈ + +L + (3.30) 34 Eq.3.30 is an example of n point approximation. Here ω1,ω2K ωn are weights and 1 2 , n ξ ξ K ξ are the Gauss points. If a one point formula is employed, the integral becomes: 1 1 1 1 I f (ξ ) dξ ω f (ξ ) − = ∫ ≈ (3.31) If a two point formula is employed, the integral becomes: 1 1 1 2 2 1 I f (ξ ) dξ ω f (ξ ) ω f (ξ ) − = ∫ ≈ + (3.32) The finite element method naturally incorporates some error as a numerical approximation. The complex continuums were modeled with many finite elements with simple shape functions to represent the element deformations (Eq.3.8). Hence, it would be undesirable to incorporate additional error in the integration of stiffness terms. It is desirable that the integration be exact. In a one point formula two parameters are considered 1 1 (ω and ξ ) . Suppose that our integration is required be exact when f (ξ ) a polynomial of is order one. So, suppose f (ξ ) is a function 0 1 f (ξ ) = a + a ξ . If we select 1 1 ω = 2and ξ = 0 Gaussian quadrature will yield an exact result. In a two point formula there are four parameters to choose 1 2 1 2 (ω ,ω ,ξ and ξ ) .Suppose that f (ξ ) must be exact for a cubic polynomial, 2 3 0 1 2 3 f (ξ ) = a + a ξ + a ξ + a ξ . The error term will be, 1 2 3 0 1 2 3 1 1 2 2 1 Error (a a ξ a ξ a ξ ) dξ [ω f (ξ ) ω f (ξ )] 0 − = ∫ + + + − + = (3.33) 35 Solution yields four nonlinear equations and they have the unique solution, 1 2 1 2 1 1/ 3 ω ω ξ ξ = = = − = − (3.34) As a result by using two Gauss points and by using the values in Eq.3.34 a cubic expression or less can be integrated exactly. By increasing the number of Gauss points, different weights (ω) and different locations (ξ ) can be found. In the FE algorithm that is developed in this study two Gauss points are used. Two Dimensional Integrals The equation 3.29 for our stiffness terms involves twodimensional integrals. So we need to extend Gaussian quadrature to the two dimensional integral form: 1 1 1 1 I f (ξ ,η) dξ dη − − = ∫ ∫ (3.35) If I is in a form like Eq.3.35, I can be written as, 1 1 1 1 1 1 1 [ ( , ) ] [ ( , ) ] ( , ) n i i i n n j i i j j i n n i j i j i j I w f d w w f w w f ξ η η ξ η ξ η − = = = = = ≈ ≈ ≈ ∫ Σ Σ Σ ΣΣ (3.36) Stiffness matrix (Eq.3.29) is twodimensional integral. The product of BTDBdet J is quadratic in terms of ξ andη .So two point Gauss Quadrature yield an exact result. It is 8x8 36 matrix so it has 64 elements. Each term must be calculated by using Eq.3.36. If ( , ) ( T det ) e ij f ξ η = t B DB J putting two for n in Eq.3.36 yields, 2 2 1 1 1 1 2 1 2 2 1 2 1 2 2 2 ( , ) ( , ) ( , ) ( , ) ij k ≈ w f ξ η + w w f ξ η + w w f ξ η + w f ξ η (3.37) where 1 2 1 1 2 2 1 1 1 3 3 w w ξ η and ξ η = = = = − = = (3.38) After input of these weights and Gauss points into the Eq.3.37 kij can be found as, 1 1 1 1 1 1 1 1 ( , ) ( , ) ( , ) ( , ) 3 3 3 3 3 3 3 3 ij k ≈ f − − + f − + f − + f (3.39) ξ η 1 1 ( , ) 3 3 1 1 ( , ) 3 3 1 1 − ( , ) 3 3 − − 1 1 ( , ) 3 3 − 1 2 3 4 Figure 3.3 Two Points Gauss Quadrature using 2x2 rule. 37 3.4 The Global Stiffness Matrix and Matrix Banding The element stiffness matrix is an 8x8 matrix and it has 64 elements. After calculating all stiffness matrices for all elements, the global stiffness matrix need to be formed from the element stiffness matrices. In Fig. 3.4 the plate is divided into n elements, this procedure is called meshing. For every element each of the stiffness terms ( , ) ( T det ) e ij f ξ η = t B DB J is evaluated using Eq.3.29. While calculating the terms of global stiffness matrix, if a node is only used by one element, its stiffness terms should be placed in the global stiffness matrix directly (like Fig.3.4 node 1). If a node is used by two elements (like Fig.3.4 nodes 3,5..) , stiffness terms for this node from two elements must be added to form the stiffness terms of the global stiffness matrix for that node. Lastly if a node is used by four elements (like Fig.3.4 nodes 4,6…) stiffness terms for that node from four elements must be added to form the global stiffness matrix. Figure 3.4 A Simple Model for Global Stiffness Matrix 38 Figure 3.5 A Quadrilateral Element and Nodal Displacements In Eq.3.40 and Eq.3.41 the stiffness matrices for element one and element two displayed in Fig.3.4 are shown 1 1 1 1 1 1 1 1 11 12 13 14 15 16 17 18 1 1 1 1 1 1 1 1 22 23 24 25 26 27 28 2 1 1 1 1 1 1 33 34 35 36 37 38 3 1 1 1 1 1 44 45 46 47 48 4 1 1 1 1 5 55 56 57 58 1 1 1 6 66 67 68 1 1 7 77 78 1 8 88 q q q q q . q q q k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k sym k k k k k k 1 2 3 4 5 6 7 8 F F F F F F F F = (3.40) 39 2 2 2 2 2 2 2 2 11 12 13 14 15 16 17 18 5 2 2 2 2 2 2 2 22 23 24 25 26 27 28 6 2 2 2 2 2 2 33 34 35 36 37 38 7 2 2 2 2 2 44 45 46 47 48 8 2 2 2 2 9 55 56 57 58 2 2 2 10 66 67 68 2 2 11 77 78 2 12 88 q q q q q . q q q k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k sym k k k k k k 5 6 7 8 9 10 11 12 F F F F F F F F = (3.41) After global stiffness matrix is formed, all displacements are calculated by using Gaussian elimination method. After finding displacements, strains and stresses for all elements can be calculated. In the algorithm that is developed the user decides the mesh density. 1000 or more elements may be needed for long spans simulations. The resulting stiffness matrix will be large and much computational time will be required for solutions. Computational times can be large because a multistep solution will be required where loads are slowly increased and the [D] matrices updated after each load step. Thus the size of the stiffness matrix becomes important because the system of updated equations will be solved many times. The stiffness matrix in our problem is a symmetric matrix. Instead of using the whole stiffness matrix the banded form of the stiffness matrix can be employed and reduce the computation time. To explain the form of the banded matrix, assume that the plate in Fig.3.4 is meshed with only two elements. In Eq.3.30 and 3.31 the stiffness matrixes of two elements were given. So, the global stiffness matrix for the plate will be: 40 1 1 1 1 1 1 1 1 11 12 13 14 15 16 17 18 0 0 0 0 1 1 1 1 1 1 1 22 23 24 25 26 27 28 0 0 0 0 1 1 1 1 1 1 33 34 35 36 37 38 0 0 0 0 1 1 1 1 1 44 45 46 47 48 0 0 0 0 1 2 1 2 1 2 1 2 2 2 2 2 55 11 56 12 57 13 58 14 15 16 17 18 1 2 1 2 1 2 66 22 67 23 68 24 k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k + + + + + + + q1 q2 q3 q4 q5 2 2 2 2 25 26 27 28 q6 1 2 1 2 2 2 2 2 q7 . 77 33 78 34 35 36 37 38 1 2 2 2 2 2 q8 88 44 45 46 47 48 2 2 2 2 q9 55 56 57 58 q10 2 2 2 66 67 68 q11 2 2 77 78 q12 2 88 k k k sym k k k k k k k k k k k k k k k k k k k k k k k k + + + F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 = (3.42) The banded form of 3.42 is shown in 3.43. The bandwidth of the matrix has been reduced to 8. 1 1 1 1 1 1 1 1 11 12 13 14 15 16 17 18 1 1 1 1 1 1 1 22 23 24 25 26 27 28 0 1 1 1 1 1 1 33 34 35 36 37 38 0 0 1 1 1 1 1 44 45 46 47 48 0 0 0 1 2 1 2 1 2 1 2 2 2 2 2 55 11 56 12 57 13 58 14 15 16 17 18 1 2 1 2 1 2 2 2 2 66 22 67 23 68 24 25 26 2 k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k + + + + + + + 2 7 28 1 2 1 2 2 2 2 2 77 33 78 34 35 36 37 38 1 2 2 2 2 2 88 44 45 46 47 48 2 2 2 2 55 56 57 58 2 2 2 66 67 68 2 2 77 78 2 88 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 k k k k k k k k k k k k k k k k k k k k k k k k k + + + (3.43) Banded matrix equation solvers exist which helps to greatly reduce the computational time associated with solving set the set of equations during each load step. 41 To conclude, the element stiffness terms are calculated from Eq.3.29 where ( , ) ( T det ) e ij f ξ η = t B DB J . After developing element stiffness matrices, the global stiffness matrix is formed. After the total potential energy is formed, which is the sum of the strain energy terms from Eq.3.28 and the work potential of internal and external forces, the theory of minimum total potential can be used. The resulting systems of equations that must be solved have the form: K q = F (3.44) This set of equations is arranged in the banded form discussed and then solved for the nodal deformations {q}. After finding displacements strains can be calculated from, ε = Bq (3.45) The stress and elastic strain components are related by a set of coefficients known as Generalized Hooke’s Law. This law can be written as, σ = Dε = D B q (3.46) Where D is equivalent elasticity matrix. 42 3.5 A [D] Matrix for Plane StressIsotropic Material and Partly Wrinkled Membrane Elements The D matrix is used while forming element stiffness matrixes and finding element stresses. The algorithm that is developed is nonlinear since the membrane D matrix is dependent on principal strains. In our models, elastic elements are used for the web on the rollers and for the exiting free span. Partly wrinkled membrane elements that were developed by Miller and Hedgepeth [5] are used for the web in test free span. The D matrix can be explicitly stated for any material. For the plane stress state for an isotropic material, stress strain relations can be written as: 1 ( ) 1 ( ) 2 (1 ) x x y y y x xy xy v E v E v E ε σ σ ε σ σ γ τ = − = − = + (3.47) By solving for the stresses in Eq.3.47 the D matrix can be found: 2 1 0 1 0 1 (1 ) 0 0 2 v E D v v v = − − (3.48) In taut membrane behavior, both principal stresses are greater than zero and the web cannot trough or buckles out of plane. In taut behavior, the two inplane principal strains may both be greater than zero or 1 ε can be greater than zero and 2 ε less than zero as long as the ratio 43 (ε 2 /ε1 ) is greater than the negative value of Poisson’s ratio ( −v ) . The [D] matrix given in Eq.3.48 is used to relate stress to strain in (1) membrane elements exhibiting taut behavior and (2) in linear elastic elements, used to model the web upon rollers. For the linear elastic elements there are no conditions placed upon the use of this [D] matrix, the principal stresses and strains can take on positive and negative values. When unsupported membrane elements have a first principal strain 1 ε less than zero this infers that the second principal strain is also negative and less than 1 ε from the rules we use to order principal strains. In this case the membrane elements exhibit a slack or unstressed behavior between stress and strain and the [D] matrix is null. In wrinkled membrane behavior the first principal stress 1 σ is greater than zero and the second principal stress 2 σ is zero. In terms of strain this behavior is represented with a first principal strain 1 ε greater than zero but now the second principal strain 2 ε is always negative. Not only 2 ε is negative, but the ratio of ( 2 1 ε /ε < −v ). For a simple membrane in uniform tension with no lateral constraint we would expect a lateral contraction governed by the expression 2 1 ε = −vε . We would expect this lateral contraction occur while the web remained planar, no buckling would be expected. However if 2 ε becomes more negative than 1 −vε we would expect a compressive 2 σ stress to have developed but since membranes by definition can withstand no compressive 2 σ stress without buckling, we would assume this element has entered the wrinkled state. For the wrinkled membrane positive nonzero principal stress can be supposed to act along the wrinkle. In a wrinkled membrane element in longitudinal tension, if 1 σ is the nonzero positive principal stress, the longitudinal direction 44 will be observed to be parallel to the wrinkles. By the property of Mohr’s circle it is known that 1 2 2 2 x y R O ε ε ε ε − = + = (3.49) where R is the radius of Mohr’s circle for strains and O is the distance from the origin of the coordinate system to the center of the circle. Consider a thin flat membrane in a state of plane stress in an xy coordinate system. Let the principal axes I and II be rotated through an angle α relative to the xy axes (Figure 3.6). 2α xy γ x ε y ε ε Figure 3.6 Mohr’s Circle for Plane Strain From the Fig.3.6 45 1 2 1 2 cos 2 2 cos 2 2 x y R R ε ε ε α ε ε ε α + = + + = − (3.50) Substituting the value of R (Eq.3.49) into Eq.3.50 yields 1 2 1 2 1 2 1 2 ( ) cos 2 2 2 ( ) cos 2 2 2 x y ε ε ε ε ε α ε ε ε ε ε α − + = + + − = − (3.51) Eq. 3.51 yields 1 2 1 2 (1 cos 2 ) (1 cos 2 ) 2 2 (1 cos 2 ) (1 cos 2 ) 2 2 x y ε ε ε α α ε ε ε α α = + + − = − + + (3.52) If P = cos 2α , then Eq.3.52 becomes 1 2 1 2 (1 ) (1 ) 2 2 (1 ) (1 ) 2 2 x y P P P P ε ε ε ε ε ε = + + − = − + + (3.53) Also xy γ can be calculated from Fig.3.6. 1 2 1 2 2 sin 2 2( ) sin 2 ( ) sin 2 2 xy xy γ R α ε ε γ α ε ε α = − = = − (3.54) If Q = sin 2α then Eq.3.54 becomes 1 2 ( ) xy γ = ε −ε Q (3.55) 46 As a result, the usual strain transformation equations yield, 1 2 1 2 1 2 (1 ) (1 ) 2 2 (1 ) (1 ) 2 2 x y xy P P P P Q Q ε ε ε ε ε ε γ ε ε = + + − = − + + = − (3.56) where 1 2 1 2 cos 2 sin 2 x y xy P Q ε ε α ε ε γ α ε ε − = = − = = − (3.57) Within a wrinkled region usual elasticity equations do not apply. Instead, the assumption of negligible bending stress in the membrane yields the stress state as: σ1 = Eε1 ; σ 2 = 0 (3.58) Or 1 1 1 1 1 1 1 1 1 1 (1 ) (1 ) 2 2 (1 ) (1 ) 2 2 1 1 2 2 x x y y xy xy P P E P P E Q QE σ σ σ ε ε ε σ σ σ ε τ σ τ ε = + = + = − ⇒ = − = = (3.59) Expressing stresses in terms of strains in the form of σ = Dε is desirable for numerical analysis. D is 3×3 matrix, ( , , ) T ( , , ) T x y xy x y xy σ = σ σ τ and ε = ε ε γ . Because the problem is statically determinant within the wrinkled region (for example 2 σ = 0 independent of the 47 valuesε1 and ε 2 ) D is singular and many possible representations for D are possible. If λ plays role of Poisson’s ratio then 2 1 λ = −ε /ε . In a wrinkled region λ is always greater Poisson’s ratio [19]. In some points λ can take the value 1. So, choosing D matrix similar to Eq.3.48 is not a right choice because of the 1/(1−λ 2 ) term. Consider D matrix is like a b c b d e c e f (3.60) If Eq.3.60 substitute into σ = Dε and impose Eq.3.56 and Eq.3.59 yields (1 ) (1 ) 2 0 0 0 (1 ) (1 ) (1 ) 2 0 0 0 0 0 (1 ) 0 (1 ) 2 0 (1 ) 0 (1 ) 0 (1 ) 2 0 0 0 0 (1 ) 0 (1 ) 2 0 0 (1 ) 0 (1 ) 2 0 P P Q a P P P Q b P P Q c P E P P Q d P P Q e Q P P Q f + − + − + − + − − = − + − + − − + − (3.61) Solutions for the elements of D matrix are not unique because the coefficient matrix in Eq.3.61 is not singular. If b=0 selected and replaced into Eq.3.61 it is found that (1 ) ; (1 ) 2 2 ; 4 4 E E a P d P E E c e Q f = + = − = = = (3.62) The resulting D matrix is the matrix which is presented by Miller and Hedgepeth [5] for partly wrinkled membranes. D matrix is: 48 2(1 ) 0 0 2(1 ) 4 0 1 P Q E D P Q Q + = − (3.63) where P and Q stated at Eq.3.57. As a result, in our models elastic elements and wrinkled membrane elements are used. For the elements in the web region on the rollers, the D matrix (Eq.3.48) is used to relate stress to strain. The elements in the web span have strain dependent D matrices as explained below. Although the slack behavior is possible in some applications it is not possible in a nonlinear formulation which employs an incremental force method. Once edge slackness begins in a web using a force method, convergence to an equilibrium solution is not possible. Convergence would be possible using an incremental displacement solution. Since this research is focused on applications where edge slackness does not occur, an incremental force solution was acceptable. 1 σ 1 σ 1 σ 1 σ 2 σ 2 σ 2 σ = 0 1 σ = 0 2 σ = 0 1 1 2 0 <ε and υε > −ε 1 1 2 0 <ε and υε < −ε Figure 3.7 Behaviors of Wrinkling Membrane Elements [12] D matrices are defined for all allowable behaviors (taut, wrinkled and slack). In Fig 3.7 behaviors of wrinkling membrane elements can be seen. A useful algorithm for choosing D 49 matrix may be expressed as [5]: 1 1 1 2 ; 0 ; 0 ; s w T D D D D and D D otherwise ε ε υε ε = < = > < − = (3.64) Where D matrices are defined as, 0 s D = (3.65) For slack behavior, as 2 1 0 1 0 1 0 0 (1 ) / 2 T E D υ υ υ υ = − − (3.66) For taut behavior, and as 2(1 ) 0 0 2(1 ) 4 1 W P Q E D P Q Q Q + = − (3.67) For wrinkled behavior where P = cos 2α ,Q = sin 2α and as stated in Eq.3.57. These algorithms will be used to establish the code that will be explained in the following chapters. 50 4. CHAPTER IV 4.1 MISALIGNED ROLLERS In this chapter the work done of Beisel [12] will be briefly reviewed. The algorithm developed will be explained step by step. Also the code that implements the algorithm will be described. Then the measures taken to decrease CPU time and to automate the code will be discussed. Finally a new slack edge criterion for misaligned rollers will be developed. 4.2 Beisel’s Method for Modeling Wrinkles Due to the Misaligned Rollers Beisel developed a method to model wrinkle formation due to a downstream misaligned roller. Like Webb he used commercial FE code COSMOS. He used wrinkling membrane elements in free spans to allow troughs to form. In Fig.4.1 Beisel’s wrinkling model for 51 misaligned roller is shown. Figure 4.1 Beisel’s FE Wrinkle Model for Misaligned Roller Contrary to Webb, he used five panels of elements rather than three panels. The first panel represents the web on the upstream roller. The second panel represents the free span where wrinkling membrane elements are used. The third panel represents the web on the downstream misaligned roller. He employed the fourth and fifth panels to enforce desired boundary conditions. First, he applied tension to the web when he reached the desired tension load he began to apply shear force to the model as shown in Fig. 4.1. Beisel employed this five panel model for the following reasons: A. The asymmetric shear forces allowed him to model the zero moment boundary condition at the misaligned roller found by Shelton [17]. B. The fourth panel was modeled using regular elastic elements. This was done because the fourth panel acts to increase the bending stiffness of the elastic elements in panel three. The failure sequence of events was: 1. Troughs form at a critical angle of misalignment given by Beisel’s previous 52 work [2]. In the finite element code the troughs are modeled by elements that may assume the wrinkled membrane behavior described in Chapter III. This happens whenever σ 2 < 0 . In a real web this does not occur until the 2 σ stress becomes more negative than the buckling stress in the free span ycr σ ( 3(1 2 ) x ycr h E a π σ σ ν = − − ; h is thickness, and a is the distance between two rollers [12]). The onset of troughs can be predicted with linear buckling theory with closed form expressions developed by Beisel [12]. To predict wrinkling of the web on rollers requires nonlinear analysis since the entering free span has already buckled in the form of troughs. 2. After troughs form CMD compressive stresses begin to appear in the elastic elements in panel three that border panel two. As the shear forces and the associated misalignment is increased further, these CMD compressive stresses become more negative and finally approach the value in Eq. 2.1, at which point wrinkles are eminent. The elastic elements in the panel four restrict the bending in panel three due to the troughs that have formed in panel two. Beisel increased the shear forces until the critical compressive stress given by Eq. 2.1 was induced in the linear elastic elements at the entry of the misaligned downstream roller. Then, he concluded that the rotation of the nodes at the entry of the downstream roller should be equal to the angle misalignment in the roller. Beisel and Webb ran similar experiments to determine the onset of wrinkle formation due to the misaligned downstream roller. Beisel compared his results with these experimental results. He modeled a polyester web with a thickness of 0.00092 in (92 gages). The web 53 parameters for this web were a Young’ Modulus of 712000 psi, a Poisson’s ratio of 0.3, a web width of 6’’ and again the thickness was 0.00092’’. The rollers had a radius of 1.45’’ and ycr σ from the Eq.2.1 was about 270 psi. He modeled 6’’; 18’’ and 30’’ span lengths and compared his results with the experimental findings. In Fig.4.2Fig.4.4 the comparison of model and experimental data is presented. 0 0.002 0.004 0.006 0.008 0.01 0.012 0 5 10 15 20 25 30 35 40 45 Tension (lbs) θcr (rad) . . experimental data FE wrinkle model Figure 4.2 L = 6’’ Span Results 0 0.005 0.01 0.015 0.02 0.025 0 5 10 15 20 25 30 35 40 45 Tension (lbs) cr (rad) . experimental data FE wrinkle model Figure 4.3 L = 18’’ Span Results 54 0 0.01 0.02 0.03 0.04 0 5 10 15 20 25 30 35 40 45 Tension (lbs) θcr (rad) . experimental data FE wrinkle model Figure 4.4 L = 30’’ Span Results As seen from these three charts, the results from FE model and experimental results show good agreement. So Beisel’s model was successful in estimating wrinkle formation due to misaligned roller for a typical web span. In his study he compared his results with Webb’s results. He claims that his model yields better results than Webb’s model. He also claimed that the assumption that Webb proposed ( cr θ is the sum of cr,trough θ and cr,wrinkle θ ) is not true. Beisel achieved good agreement with experimental results without relying upon Webb’s assumption. 4.3 A New Algorithm for Predicting Wrinkles Due to the Misaligned Rollers As mentioned in the introduction, there are three types of wrinkles: MD wrinkles, CMD wrinkles and shear wrinkles. Shear wrinkles can occur due to roller imperfections such as 55 misaligned rollers or tapered rollers. The goal of this new algorithm is to codify and automate the analysis that Beisel perfected using a commercial finite element code (COSMOS). The term “shear wrinkle” resulted from the realization that these troughs and finally wrinkles were the result of shear forces in the web. Both the misaligned roller and the tapered roller induce shear in the web. Beisel [12] was successful to develop a method for predicting web wrinkles on rollers by using membrane elements described by Miller and Hedgepeth [5]. He applied this method to the prediction of wrinkles due to misalignment in rollers, tapered rollers and crowned rollers and he confirmed his results with laboratory tests. Beisel used commercial finite element code COSMOS to apply this method. In this method while the elements representing the web on rollers are modeled with elastic elements, the web in the free spans are modeled with wrinkle membrane elements. These elements cannot react compressive stresses and they can be in one of three states. These states include taut web, wrinkled web and slack web. In the taut web state, the elements can resist tensile stresses in both principal directions. In the wrinkled web state, membrane elements can withstand tensile principal stress in one direction and zero stress in the other principal direction. In the slack web state, the elements can carry no stresses in any direction. In this algorithm, forces are increased in to the model step by step. In the first step, all elements are modeled with elastic elements. After the first step, the principal strains for each element in free spans are calculated and stored. By using principal strains, code will select which of D matrices given below for three states will be used for the next load step. 56 1 1 1 2 ; 0 ; 0 ; s w T D D D D and D D otherwise ε ε υε ε = < = > < − = (4.1) where D matrices are defined at Eq.3.653.67. After convergence is obtained in each step, the compressive CMD stresses in the linear elastic elements are reviewed. If those stresses remain greater than the Timoshenko shell buckling stress (Eq.2.1) the shear force would be increased. If the CMD stresses in these elements became more negative than the Timoshenko shell buckling stress (Eq.2.1), the shear force would need to decrease and a bracketing method would be employed to determine what shear force would produce a negative a negative CMD stress essentially equal to the Timoshenko buckling stress. Once accomplished, the misalignment or taper that induced that level of shear force would be determined. By using the method explained above, a finite element code will be developed in Excel VBA (Visual Basic Excel). This code can be executed in any PC with Excel installed without need of a commercial FE code license and it allows users to analyze the misaligned roller case using a simple Excel based interface. The advantage of this code will be that users will not need any linear or nonlinear finite element background to execute the code. They do not have to know the kinetic and kinematic boundary conditions for misaligned rollers. The inputs will include parameters such as web tension, web width, span length, roller diameter, Poisson’s ratio and elastic modulus. When executed the code will automatically form a finite element mesh based upon the inputs with elastic quadrilateral elements representing the web supported by rollers and with wrinkle membrane quadrilateral elements representing the web 57 in the free span. The first code will implement boundary conditions for a web approaching a misaligned roller. Other boundary conditions will be studied later. Beisel [2, 12] and Webb [9] studied the misaligned roller case. The boundary conditions that they used are first proposed by Shelton [14] and then Good et al. [4]. They considered the web span as a beam. A classic beam is one which the web span length would be ten times longer than the width. Shear effects become “important” when L 10 w < . Tension becomes “important” when the lateral deformations become large. “important” in this context means that these effects have sizable influence on the lateral deformations of the web. The boundary conditions that are used by Beisel [2, 12], Webb [9], and others will be used in this model. The validity of using these boundary conditions was verified by comparison to experimental results by these authors. The normal entry condition of a web approaching a roller was enforced using coupling equations which enforce multipoint constraints. Lines of adjacent nodes crossing a roller in panel one and panel five had their CMD displacements coupled. Each adjacent line of nodes was coupled separately and in this way Poisson contraction of the web could occur unimpeded. There was no coupling of nodes in the web in contact with the misaligned roller. Since the moment in the web in the vicinity of the misaligned roller is small or zero the deformations of nodes are nearly that associated with a rigid body rotation. This results in the normal entry condition being satisfied in the web at the misaligned roller without enforcing a 58 multipoint constraint. The lines of nodes in the CMD at the exit of panel one and entry of panel five were each coupled in machine direction displacement. This procedure was done to ensure the maximum moment in the free web spans occurred at the border between span one and span two and the border between span four and span five. The system that is modeled is shown in Fig .4.5. Figure 4.5 The System That is Modeled The system of five panels is shown in Fig.4.5 modeled in Fig.4.6. The coupling discussed earlier is also shown. 59 Figure 4.6 Misaligned Roller FE Wrinkle Model The model is divided into five sections. The first panel represents the web on the upstream roller. The coupled nodes in this panel are used to enforce normal entry and exit on upstream roller. The second panel represents the entry web span to the misaligned roller. Here wrinkling membrane elements are used to simulate web behavior which allow troughs. Different from Beisel’s model at the first attempt the fourth panel is also modeled with wrinkling membrane elements. Shear forces are applied to the web on the upstream and downstream rollers to simulate the shear, moment, and lateral deformations of a web passing over a misaligned roller. The third panel represents the web on the downstream misaligned roller. A central node is fixed in the MD and CMD directions to prevent rigid body translations of the model. Rigid body rotation is prevented by the coupling of the CMD deformations of the lines of nodes crossing panel one and panel five. The fourth panel and the fifth panel elements and boundary conditions help ensure the zero moment boundary condition at the misaligned roller. The flow chart for the program is shown in Fig. 4.7. 60 61 Figure 4.7 Flow Chart for New Algorithm 62 The stiffness terms for the elements on the upstream roller are calculated by using DT. Then these stiffness terms are assembled into the global stiffness matrix. The same procedure is followed for the other web regions on the rollers. For the web spans D matrix for the elements will be formed differently. At the first load level and first iteration DT is used for all elements. Then, for all remaining load levels and for all iterations, the program selects one of the three D matrices from (Eq. 3.64) Dtaut, Dwrinkled or Dslack by evaluating the principal strains calculated in a previous step as explained in Eq. 3.6567. The nonlinearity for this case is due to the variable D matrices for the span elements as the shear loads increase. After selecting D matrices, elemental stiffness matrices are formed and the stiffness terms will be assembled into the global stiffness matrix using the same procedure as the web on the rollers. After the global stiffness matrix is formed, lines of nodes in the MD on the first and fifth panels are coupled in the y direction and the point at the center of the model is fixed in x and y directions as shown in figure 4.6. Next, the shear and traction forces are applied to the system. From the set of equations KQ=F the displacements {Q} can now be calculated. Strains {ε } and stresses {σ } can then be determined using the displacements. The strains are calculated usingε = BQ , stresses are calculated with the aid of σ = DBQ for all elements. After calculating strains, principal strains are also determined from the cartesian strains so that the proper D matrices can be selected for the next iteration or for the next load level. The shear and tension forces are applied incrementally in 20 steps. Five iterations are performed for each load step to allow P and Q to converge. As mentioned before if the 63 principal axes one and two are rotated through an angle α relative to the xy axes, P is the cosine of that angle and Q is the sine of that angle. For each new load step, the code first uses the principal strains calculated in the 5th iteration of the previous load step to determine the state of the wrinkling membrane for the 1st iteration. Then, the strains calculated in the first iteration for that load level are used to select D matrices for the second iteration. The procedure will continue in the same way. At each iteration, the strains calculated for the previous iteration are used to select D matrices till the maximum iteration number (5) is reached for that load step. The analysis proceeds in load steps with iterative analysis steps occurring within each load step. The iterative analysis steps are necessary to allow the values of P and Q to converge in the elements with the wrinkled behavior prior to moving to the next load step. In Fig.4.8 it shows how P and Q (Eq.3.57) behave with iteration for an element in Panel 2. Figure 4.8 Convergences of P and Q 64 Figure 4.9 Excel Input for the Program In Fig. 4.9, the screen where the web parameters are input to the code is shown. As seen from the Excel input screen, the user enters the parameters such as width of the web, a roll dimension (quarter circumference of the roller), a span length dimension, the thickness of the web, the elastic modulus of the web, Poisson’s Ratio of the web, and the web tension in units of traction (stress) in the x direction. Also the shear in y direction in units of traction (stress) to the Excel sheet is input. In the “m height elements” cell, user enters the mesh density along x direction. In “n1 roll elements” cell, user inputs the mesh density of the rollers along y direction. Similarly in n2 span elements cell user decides about the mesh density of the span elements along y direction. This is not the final form of input. In the next chapters the efforts for the final form of input will be discussed. 65 The output of the code developed to predict wrinkling due to downstream misaligned roller are σ x and σ y stresses for all elements on upstream roller, the downstream roller and in the web span between the upstream and the downstream roller are shown in Figure 4.10. Figure 4.10 Excel Output for Stresses For this case critical buckling stress by Eq.2.1 is about 270 psi. The marked row of stresses are the stresses of the elements at the entry of the downstream roller. The elements at the entry of the roller buckle first. After the output is displayed the minimum (most negative) stress in these elements should be compared with the critical shell buckling stress. If unequal, then the shear force should be increased or decreased until the minimum compressive stress 66 in the elements at the entry of the downstream roller reaches the critical shell buckling stress. In Fig.4.10, the maximum compressive stress (269.6) reached the Timoshenko shell buckling stress (270). At this instant, the angle displayed as an output in Figure 4.9. is the critical angle of misalignment of the downstream misaligned roller for the onset of wrinkling (0.0073 rad.). 4.4 Comparison of Results with Previous Works The code that was developed was executed for some cases defined by Beisel and the code results are compared with Beisel’s experimental and FE results. As mentioned in the previous chapter he used a commercial FE code called COSMOS and modeled the misaligned roller case. He performed his tests using 92 gage polyester web to verify his model. For this web, Young’s Modulus is 712000 psi, Poisson’s ratio is 0.3, the width of web is 6’’ and the thickness is 0.00092’’. The roller has a radius of 1.45’’. Eq.2.1 yields about 270 psi for the critical shell buckling stress for this case. Beisel was comparing the stresses at the nodes at the entry of the roller with the critical compressive stress predicted by Timoshenko shell buckling criteria. In this study, compressive stresses in the elements at the entry of the roller are used in the comparison instead of nodal stresses. Stresses at four Gauss points of the elements are computed and their averages are found as elemental stresses. The following graphs in Fig. 4.1113 show the experimental critical angle of misalignment at 67 the downstream roller, the angle predicted by Beisel’s model and by Yurtcu’s model for 6’’, 18” and 30’’ web span lengths. 0 0.005 0.01 0.015 0.02 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester,L=6" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu's Model(VBA) Figure 4.11 Comparison of 6’’ Span Results 0 0.005 0.01 0.015 0.02 0.025 0.03 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester, L=18" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu's Model (VBA) Figure 4.12 Comparison of 18’’ Span Results 68 0 0.01 0.02 0.03 0.04 0.05 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester,L=30" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu's Model (VBA) Figure 4.13 Comparison of 30’’ Span Results From the figures, it can be seen that the results of the code developed agree well with the lab test data and the results of the commercial code. 4.5 Improvement of Execution Time It was found that execution the code for long spans that twenty load steps with five internal iterations in each load step required a large amount of CPU time. For example, 2000 elements may be employed in a long span case. For every element we have an (8*8) stiffness matrix. Forming the global stiffness matrix for 2000 elements and solving it 100 times (five iterations in 20 load steps) required extensive CPU time. Although we use banded matrix in 69 our code long spans could require three hours to produce a result on a single core computer. Due to these long execution times the next focus was decreasing the CPU time. First the number of iterations was decreased from five to two, then from two to one, and the result was not significantly changed. Then step by step the load step increments were decreased from twenty to four, and still reasonable results were achieved. However, it was found that decreasing the load step increments to less than four caused the results to err dramatically. This problem is appears to converge with four load increments and one iteration within each load increment. It was mentioned in the previous section that Beisel modeled the second span with elastic elements, but that wrinkling membrane elements were used in this code. This means that we were calculating principal strains to select which D matrix would be used for every element in span two. Beisel ran taut elements in the second free span because the model was a 900 wrap case. For the misaligned roller case this subjects the upstream span to shear and the downstream span to twist. A web span will absorb a large amount of twist without forming negative 2 σ stresses. So Beisel assumed the D in all elements in the downstream span would remain taut. We forced D matrices to be DT in the second free span. This also provided better agreement between our code and other results. Changing element types from wrinkling membrane elements to elastic elements also helped to decrease CPU time. Our code produced reasonable results, and did this within seconds. After making these changes, the flowchart of our program will be like Fig.4.14. The name “WRINKLINGsystem” will be applied to the part of the new chart that begins after “Mesh Model with Quadrilateral Elements” and continues to “Load Level<4”. This name will be used while attempting to automate the code in Part 4.6. 70 71 Figure 4.14 Flow Chart of Computer Code for Improved Execution Speed 72 4.6 Automating the Code From Fig 4.9 it can be seen that the user enters parameters such as web width, roller wrap arc length, span length, web thickness, web elastic modulus, Poisson’s ratio, the number of elements across the web width, the number of elements down the span length, the number of elements across the roller, the web tension and shear. The research objective was to limit the inputs that could be comprehended by users with no knowledge of the finite element method. These would include only web width, roller wrap arc length, span length, web thickness, web elastic modulus, Poisson’s ratio and web tension. One of the inputs required by the code is web tension. Web tension is one of a very few parameters that are controllable in a web process machine. Thus it would be optimal for a chart to be produced for the user that shows how much misalignment is allowable as a function of web tension, rather than computing what misalignment in a roller is acceptable at one tension. The user can then decide to solve an instability problem by better aligning the rollers or by changing the web tension. Thus the user should enter only certain parameters and the code will determine the mesh parameters and shear force required to induce wrinkling. To automate the mesh parameters the code was executed several times to explore what mesh density was required to produce a converged result. Very reasonable results were obtained by dividing the width and span length per piece per dimension and dividing the rollers with the integer part of six times of onefourth of roller circumference. The mesh density is delimited for very long span lengths, very short span lengths, very long wide webs and very short wide webs. The meshing procedure and convergence check will be addressed more detailed at 73 Chapter VIII. The next step was determining how to automate the search for the shear force that was required to induce the wrinkle instability. The user inputs the web tension. A linear interpolation scheme was used to determine the shear force. One level of shear will produce a certain level of compressive y σ stress at the misaligned roller, which can be determined by the code. A second level of shear will produce another level of compressive y σ stress. Interpolation can be used to estimate the level of shear that will produce the shell buckling stress. It is an estimate because this is nonlinear analysis. That estimate can then be input to the code to help refine the actual shear level that will induce wrinkling. A slack edge criteria was used as a starting point. If a slack edge forms during the computations the code will fail. This is because increased shear will not result in increased compressive stress in the web on the roller, it will only increase slackness. For more information about slack edge criteria earlier work done by Good [40] can be visited. Web spans can be modeled as beams [12]. From EulerBernoulli beam bending theory the stress in the crosssection is: F My A I σ = ± (4.2) If this stress is equal to zero, a slack edge occurs. If our web has a thickness of t, width of w, span length of L and the applied traction in x direction is Tx the shear stress (Txy) for the web to be slack can be calculated from Eq.4.2 as: 0 2 x w M T I σ = = − (4.3) Where moment M is calculated from: 74 M = Txy ×w×t × L (4.4) If this value put into Eq.4.2 and do the calculations Txy (Tslack) will be found as 6 x xy T w T L × = (4.5) This was used as the starting value of Traction XY (shear) in the code (Fig.4.9).It found that for all cases computed the compressive stress in the first row of elements in the roller two was lower (less negative) than the value that was calculated from Timoshenko buckling criteria (Eq.2.1). Thus it found that if the shear value in Eq.4.5 was used as traction xy (Tslack) and half of it ( Tslack/2) one could be sure that there were two data points which would produce compressive stresses less than that is needed to buckle the web (Tcritical) (Figure 4.15). critical σ 2 σ 1 σ Figure 4.15 Linear Interpolation I 75 After finding two points we can estimate Tcritical from the tangent tanα : tan 2 1 2 ( / 2) critical Slack Slack Critical Slack T T T T σ σ σ σ α − − = = − − (4.6) from here Tcritical can be found as: 2 2 1 ( ) ( ( / 2)) critical Slack Slack ( / 2) Critical Slack T T T T σ σ σ σ − × − = + − (4.7) The algorithm for Linear Interpretation I can be seen in Figure 4.16. The code takes the value of web tension (traction x) from Excel input page and calculates F1, F2 and F3 values with the help of Tslack (Eq.4.5). Here F1 is equal to Tslack /2, F2 is equal to Tslack and F3 is equal to 2* Tslack .Than code runs for F1 and F2 values within the WRINKLINGsystem and calculates Sigma1 and Sigma2 values. By using Eq.4.7 code calculates Tcritical . The code checks whether Tcritical is larger than two times of Tslack . This was done because during the runs it was observed that if this was not done, the value of Tcritical increases dramatically because of the angle between two points. Taking this step helps control Tcritical value. Then the code runs WRINKLINGsystem for Tcritical value and finds SigmaS. The code continues this process until SigmaS is bigger than the Timoshenko buckling criteria (Sigmacritical). At the end of Linear Interpolation I, there is a F1 value which is less than the traction that is needed to buckle the web, and there is a F2 value which is more than the traction that is needed to buckle the web (Figure 4.17). 76 Figure 4.16 Flow Chart for Linear Interpolation I 77 The code for linear interpolation one will be: TRACTIONXY = (TRACTIONX * H) / (6 * (L1 + L2 + L1 + L2 + L1)) F2 = TRACTIONXY F3 = 2 * TRACTIONXY Do Until sigmaS > syrc FORCEXY = F2 * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigma2 = critical F1 = F2 / 2 FORCEXY = F1 * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigma1 = critical FS = ((sycr  sigma2) * (F2  F1)) / (sigma2  sigma1) + F2 If FS > F3 And control = 0 Then FS = F3 control = control + 1 End If FORCEXY = FS * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigmaS = critical F1 = F2 F2 = FS Loop 78 After establishing values F1 and F2 we can start Linear Interpolation II process. The new problem will look like Figure 4.17. critical σ 2 σ 1 σ Figure 4.17 Linear Interpolation II From Fig.4.17 Tcritical can be found as: 1 2 1 ( ) ( 2 1) critical 1 critical F F T F σ σ σ σ − × − = + − (4.8) Here the aim is to approach Tcritical value by changing the values of F1 and F2. The flowchart for Linear Interpolation II can be seen at Figure 4.18. 79 Figure 4.18 Flow Chart for Linear Interpolation II From Linear Interpolation I the values F1 and F2 are known. By using Eq.4.8 Tcritical was calculated and SigmaS was found. If SigmaS is greater than Sigmacritical we replace Tcritical with F2. If not we replace Tcritical with F1. This process was continued until SigmaS is 80 between the limits. In case of misaligned roller case the lower limit was set to 0.99*Sigmacritical and the upper limit was set to 1.01*Sigmacritical. The code for Linear Interpolation II will then be: Downlimit = 0.99 * sycr Uplimit = 1.01 * sycr Do Until Downlimit <= sigmaS And sigmaS <= Uplimit FORCEXY = F2 * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigma2 = critical FORCEXY = F1 * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigma1 = critical FS = ((sycr  sigma1) * (F2  F1)) / (sigma2  sigma1) + F1 FORCEXY = FS * H * te NFORCEXY = FORCEXY / m Call WRINKLINGsystem sigmaS = critical If sigmaS > sycr Then F2 = FS ElseIf sigmaS < sycr Then F1 = FS End If Loop Linear Interpolation I, Linear interpolation II and WRINKLINGsystem can be seen in the 81 Appendix. After automating, Excel input of the code will look like Figure 4.19. Figure 4.19 Excel Input for the Automated Program Here the user is supposed to enter web width, span length, thickness, elastic modulus, Poisson’s ratio, roller radius and the web tension. Other parameters are calculated automatically. Meshing elements are calculated with Excel equations and shear traction xy required to induce wrinkles is calculated within the code as mentioned above. After execution the code provides the following output: total time of execution (67 seconds), the maximum compressive stress (266.8 psi) in the first row of elements and the roller misalignment angle (in degrees) that produced that compressive stress. For this case the sigma critical value calculated from Timoshenko buckling criteria (Eq.2.1) is 270 psi. So the interpolation scheme discussed has produced a compressive stress in the first row of elements 82 that is very close to the critical value. The results from the automated code are shown in Figures 4.2022 and the results are compared with Beisel’s commercial FE results and his test data. 0 0.005 0.01 0.015 0.02 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester,L=6" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu`s Modified Code (VBA) Figure 4.20 Comparison of 6’’ Span Case with the Modified Code 0 0.005 0.01 0.015 0.02 0.025 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester, L=18" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu`s Modified Code (VBA) Figure 4.21 Comparison of 18’’ Span Case with the Modified Code 83 0 0.01 0.02 0.03 0.04 0.05 0 2000 4000 6000 8000 Roller Misalignment (Rad) MD Tension (psi) Wrinkles Due to the Misalignment,92 ga Polyester,L=30" Beisel`s Model (COSMOS) Beisel Lab. Test Yurtcu`s Modified Code (VBA) Figure 4.22 Comparison of 30’’ Span Case with the Modified Code The execution time of the modified code is much less than the previous version. In the previous version execution time was around two to three hours and code was able to give one result for the specified shear force. For the modified code the execution time is around five to ten minutes and the modified code finds the right shear force that will buckle the web by itself. Over the parameter ranges of these examples the modified code appears to mesh the problem adequately and yields good results at all tension levels. 84 4.7 A New Slack Edge Criteria for a Misaligned Roller In Figure 4.23, buckling region of a 92 gage polyester web can be seen. For this web Young’s Modulus is 725000 psi, Poisson’s ratio is 0.3, the width of web is 6’’, length of the web is 40” and the thickness is 0.00092’’. If 6000 psi MD stress applied to the system we are able to see wrinkles after the trough formation. If 2000 psi applied to the system after the formation of troughs slack edge occurs and we are not able to see wrinkles. Figure 4.23 Slack Edge The code was delivered to the Web Handling Research Center sponsors. One of the sponsors used the code with a low span (L/W) ratio. In that case W (width of the web) was five times 85 larger than L (length of the web). With that specific web tension that was applied to the web a slack edge was supposed to occur. But Excel VBA code was came up with a result which means wrinkles was occurred after the trough formation instead of a slack edge. The case was also modeled with COSMOS and the result was very similar with Excel VBA code. In the derivation of Eq.4.5 shear deformation was neglected. For L/W values of 0.2 it was obvious that shear deformation is important. Thus it was determined that a new slack edge criteria was needed that did account for shear deformation. Figure 4.24 Slack Edge Criteria The previous slack edge expression was derived based on Euler beam theory: ( ) x SlackEdge T a b Ehb θ = (4.9) 86 Here, a is length of the web, b is width of the web (Fig.4.23), h is thickness and E represents elastic modulus. To derive a new criteria the Timoshenko beam theory was used that incorporates both shear and tension stiffening. Figure 4.25 Beams with Shear If w is deformation, the slope of the beam isα , and can be found from: ' w wb ws w x x x α θ γ ∂ ∂ ∂ = = = + = + ∂ ∂ ∂ (4.10) Here θ is related to moment and γ is related to shear. ( ) ( ) s M x F x and x EI GA θ γ ∂ = = ∂ (4.11) The moment can be calculated from Eq.4.10 by using Eq.4.11 ( ) 1 " , , . x x s M x F Moment w EI GA x θ γ ∂ = = + = + ∂ (4.12) 87 Figure 4.26 Moment in a Beam At the boundary x = 0 , 0 0 x x S w F x GA γ = = ∂ = = ∂ (4.13) The moment M(x) = F(L − x) and F(x) = F , and substituting into Eq.4.12 yields ( ) " 0 F L x w EI − = + (4.14) Integrating yields 2 1 ' ( ) 2 F x w Lx C EI = − + (4.15) From Eq.4.13 1 x 0 S F C GA γ = = = so Eq.4.15 becomes 2 ' ( ) 2 S F x F w Lx EI GA = − + (4.16) Finally integration of Eq.4.16 yields, 88 3 2 2 ( ) 6 S F x Fx w Lx C EI GA = − + + (4.17) It is known that 0 0 x w = = . So 2 C becomes 0. We seek critical θ at x=L .From Eq.4.16, 2 1 ' ( ) x L critical 2 S L w F EI GA θ = = = + (4.18) Solving for the steering force F yields: 2 1 ( ) 2 critical S F L EI GA θ = + (4.19) This allows us to calculate the moment at x=0, 2 . 1 ( ) 2 critical S M F L L L EI GA θ = = + (4.20) The beam bending stress from moment can now be calculated as, max 2 1 1 2 ( ) 2 B critical S M w y L I L I EI GA θ σ = = + (4.21) The normal stress due to web tension can be calculated from, T T wb σ = (4.22) In taut behavior σ T >σ B if σ T =σ B the slack case occurs. So if σ T =σ B , 89 2 1 2 ( ) 2 critical S L w T L I wb EI GA θ = + (4.23) From here critical θ can be calculated as, 2 2 2 ( ) critical S T L I Lw b E GA θ = + (4.24) For a rectangular cross section 3 5 , 12 2(1 ) 6 S bw E I G and A bw v = = = + . Substituting into Eq.4.24 yields, 2 2 2 (1 (1 )) ( ) 5 critical TL w v w Ebw L θ = + + (4.25) In the calculations beginning from Eq4.10, L for span length, w for span width and t for thickness were used. If we replace them with a for length, b for width and h for thickness, it will be found, 2 2 (1 ( ) (1 )) ( ) 5 critical Ta b v b Ehb a θ = + + (4.26) If this result is compared with Eq. 4.9 it will be found that 2 2 (( ) (1 )) ( ) 5 Ta b v b Ehb a + (4.27) is the effect of shear force with Timoshenko beam theory. If v ≈ 0.3 then the term is 90 2 1 (1 ) 5 2 + v ≈ . So Eq.4.26 can be simplified as, 2 1 (1 ( ) ) ( ) 2 critical Ta b b Ehb a θ = + (4.28) If W >L (b>a) then shear effects are becoming significant. For instance if W=5L as it was in the case that we came up with earlier then, 1 (1 25) ( ) 2 13 ( ) criticalTimoshenko criticalEuler Ta b Ehb Ta b Ehb θ θ + ≅ ≅ (4.29) Eq. 4.26 supersedes Eq.4.9 and can be used for both long and short spans. 91 5. CHAPTER V 5.1 TAPERED ROLLERS In this chapter, a brief definition of tapered rollers, previous work modeling with tapered roller cases with commercial finite element codes will be reviewed. Then a new model based on using Excel VBA will be given. 5.2 Description of Tapered Rollers and Beisel’s Method for Modeling Wrinkles Due to the Tapered Rollers A tapered roller is defined as a roller with a linearly varying radius across its width [12]. Tapered rollers are commonly seen in the web handling industry. The process of roller manufacture will almost certainly result in rollers with a slight taper. These tapered rollers 92 are produced unintentionally; the use of tapered rollers in web lines may result in web damage. To solve this problem, machining techniques that involve feedback can be used, but this will be costly. Therefore, knowing the amount of taper that will not result in harm to the web process is beneficial for the industry. In Figure 5.1 a web approaching the tapered roller is shown. In the figure, the taper is somewhat overemphasized so that it can be seen. y r Roller radius Position across the web Average radius of roller Ro Figure 5.1 Tapered Roller Profile The radius of the roller at any point across the width is: r( y) = my + Ro (5.1) Here m is the slope of the roller. The velocity across the roller width is: ( ) ( ) ( ) o V y = r y ω = my + R ω (5.2) Here ω is the angular velocity of the roller. The average web velocity can be found by using average roller radius as: 93 Vavg = Roω (5.3) Variation in the velocity across the web width will cause a stress and strain upon the web: ( ) ( ) avg ( ) ( ) md md avg o o V y V my Emy y and y E y V R R ε σ ε − = = = = (5.4) These equations assume the web and roller achieve the same velocity at contact. The variation in stress across the web width causes a steering moment on the web: ( ) 2 2 2 3 2 2 12 b b j b b o o Emhy mEhb M y hydy dy R R σ − − − − = ∫ − = ∫ = (5.5) where m is roller taper, R0 is roller radius, E is Young’s Modulus, h is thickness and b is web width. Beisel [12] determined that wrinkle formation due to a downstream tapered roller is similar to the wrinkle formation due to a misaligned roller. As in the misaligned roller case, he used Timoshenko buckling criteria (Eq.2.1) to decide whether the web on the tapered roller will wrinkle or not. His wrinkle model for the tapered roller is shown in the Figure 5.2. Figure 5.2 Beisel Tapered Roller Model [12] Web span being studied Wrinkling membrane elements Nodes coupled along dotted lines, to provide normal entry, exit and travel across rollers Horizontal dotted lines coupled in y direction Vertical dotted lines coupled in x direction Uniform shear force (2fyj) applied at both ends of middle roller and end of web spans Web line tension applied at both ends of web Nodes coupled in y direction to 94 As in the misaligned case, he divided the model into five panels. The first panel is the web on the upstream roller. The second panel represents the web span being studied. The third panel represents the tapered roller. The fourth panel represents the web between the tapered roller and third roller, and the fifth panel represents the third roller. He locked the nodes in the first panel along horizontal lines. Each node on the line must move with the same y displacement as the rest of the nodes on the same line. This allows for Poisson contraction due to web line tension but requires a zero slope at the beginning of the troughed web span. He also locked the nodes at the exit of the upstream roller in the x direction along a vertical line to simulate the roller gripping the web in a noslip condition. He used wrinkling membrane elements in the second panel. Along the right edge of the second panel he locked the nodes along horizontal lines in the y direction for a very short distance. He did this to ensure the normal entry of the web to the tapered roller. He modeled the right hand side of the model to enforce the boundary conditions. He executed the model by first applying web line tension and then increasing shear force until the compressive stress across the first row of elastic nodes on the elastic wrinkling membrane element boundary reached the critical value predicted by Eq. 2.1. Then he calculated the moment associated with the last row of span elements by using Eq. 5.5. He calculated the critical taper that would induce wrinkles. He obtained the experimental results for the onset of wrinkles due to a downstream tapered roller. He compared his experimental results and model results for two materials. In the following chapters, we will compare our model results with his experimental results and with his model results. 95 5.3 A New Algorithm for Predicting Wrinkles Due to the Tapered Rollers As mentioned earlier, shear wrinkles can occur due to tapered rollers. Beisel [12] successfully modeled the tapered roller case by using the commercial finite element code COSMOS. As in the misaligned case he used wrinkling membrane elements while modeling the elements between the upstream roller and the tapered roller. The goal of our new algorithm related to the tapered roller is to codify and automate the analysis that Beisel did with the commercial finite element code COSMOS. Similar to the misaligned roller case a finite element code that calculates critical taper will be developed in Excel VBA by using wrinkling membrane elements. The boundary conditions will be similar to Beisel’s model boundary conditions. The problem is modeled material on the upstream roller, the upstream span, the web on the tapered roller, the downstream span, and finally material on the downstream roller. The system that is modeled is shown in Fig.5.3. 96 Figure 5.3 The Figure for Tapered Roller The system with a tapered roller shown in Figure 5.3 is modeled like the Figure 5.4 shown below. Figure 5.4 Tapered Roller FE Wrinkle Model 97 The new model is similar to the misaligned roller case. The model was divided into five panels. The first panel was the panel on the upstream roller. The second panel was representing the web between upstream roller and tapered roller. Here wrinkling membrane elements were used; in the rest of the model elastic elements were used. The third panel represents the tapered roller. The fourth and fifth panels help to achieve the no moment boundary condition at the tapered roller. Boundary conditions and loads were enforced which Beisel found to be appropriate for the tapered roller. The center is pinned to prevent rigid body motions. Multipoint constraints were applied on each row of nodes on the entering and exiting rollers. Along the right edge of the second panel (span one), Beisel locked the nodes along horizontal lines in the y direction for a very short distance to ensure normal entry. In contrast to his method, high mesh was used at the last row of the span one element (ten elements were used in the last row elements) and six points along the horizontal lines were locked in the y direction. This was decided after trying many ways to achieve normal entry to the tapered roller. The flowchart for the program is shown in Fig.5.5. It is similar to the misaligned roller flowchart so it will not explain it here in detail. 98 99 Figure 5.5 The Flow Chart for Tapered Roller 100 In the tapered roller flowchart the forces are applied in five time steps. After five steps from the first row of elements the moment is calculated and by using Eq. 5.5 critical taper is also calculated. Assume that the stresses of the first row of elements at tapered roller look like Fig.5.6. Here, one half of the width is demonstrated. hi Ai n ı i Figure 5.6 Calculating Moment and Critical Taper Moment is M = F h and / x σ = F A. From these two equations the moment of an element can be calculated from 1 n i i i i M σ A h = =Σ . If all moments are added the total moment of the first row of the elements of the tapered roller is found. If this value put into Eq .5.5 the critical taper for that specific case is found. At the end of the code, the code calculates total moment and taper for that case. The name “SystemWrinkler” will be applied to the part of the flow chart (Fig.5.5) that begins after “Mesh Model with Quadrilateral Elements” and continues to “Load Level <5” .This term will be used in the following chapter. 101 5.4 Automating the Code Similar to the misaligned roller case, it was aimed for the user to enter only certain parameters such as web width, roller wrap dimension, span length, web thickness, elastic modulus, Poisson’s ratio, and web tension, and run the code. The attempts to automate the code began with trying to automate the mesh. After running many cases, it was determined that reasonable results were achieved by using three elements per dimension for the web width, one element for every dimension of the web spans and the integer part of four times of the roll dimension for rollers. Excel equations were used to set these values. The meshing procedure and convergence check will be addressed more detailed at Chapter VIII. The flowcharts for Linear Interpolations I and II are used for automating the tapered roller case. Similar methods to the misaligned case were applied. Tslack and Tslack/2 (Eq.4.5) were used as a starting point for Linear Interpolation I. In Fig.4.16 and Fig.4.18, if we use the term SystemWrinkler instead of WRINKLINGsystem, the way the new flow chart works can be explained. After automating our tapered roller code the Excel input of the code will look like Figure 5.7. 102 Figure 5.7 Excel Input of Tapered Roller Here the user is supposed to enter web width, span dimension, thickness, elastic modulus, Poisson’s ratio, roller radius and the web tension. After entering these parameters with the help of the Excel equations the mesh parameters, one fourth of the circumferences of the rollers (L1 roll dimension in Fig.5.7) and the stress calculated from Timoshenko’s buckling criteria (Sigma Critical in Fig.5.7 that is found from Eq. 2.1) can be calculated. After clicking the EXECUTE button, the code runs Linear Interpolation I and Linear Interpolation II. As a result of using the code, the maximum stress at the first row of elements (MAX sigma y), total moment at the first row of elements (Moment), and the critical taper (mcr) that will result in a wrinkle for that specific element for that specific case are given as an output. 103 Figure 5.7 shows that for that element the stress calculated from Timoshenko’s buckling criteria is 247 psi. After linear interpolations, the code finds a maximum stress of 245 psi, which is very close to the Timoshenko buckling criteria. The critical taper for the case shown in Fig. 5.7 is found to be 0.00127 (in/in). The code calculated this case within six minutes. 5.5 Comparison of Results with Previous Works for the Tapered Roller Beisel [12] performed experiments and obtained data and compared his experimental results with his model. The results from the Excel VBA code will be compared with his FE model results and his experimental results. The first web he tested was a 92 gage (0.00092”) opaque polyester with a Young’s Modulus of Ex = 712000 psi, a Poisson’s Ratio of ν = 0.3, and a width of W = 6”. The nominal radius of his tapered roller was Ro = 1.49”. The compressive stress required by Eq.2.1 was approximately 265 psi for this case. Comparison of our results with his FE model and his experimental results are given below. 104 0 0.0005 0.001 0.0015 0.002 0 10 20 30 Mcr (in/in) Tension 92 ga Polyester, L=10" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.8 Wrinkles Due to Taper, 92 ga Polyester, L=10” 0 0.0005 0.001 0.0015 0.002 0 5 10 15 20 25 30 Mcr(in/in) Tension 92 ga Polyester, L=20" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.9 Wrinkles Due to Taper, 92 ga Polyester, L=20” 105 0 0.0005 0.001 0.0015 0.002 0.0025 0 5 10 15 20 25 Mcr(in/in) Tension 92 ga Polyester, L=30" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.10 Wrinkles Due to Taper, 92 ga Polyester, L=30” 0 0.0005 0.001 0.0015 0.002 0 5 10 15 20 25 Mcr(in/in) Tension 92 ga Polyester, L=40" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.11 Wrinkles Due to Taper, 92 ga Polyester, L=40” 106 Figure 5.811 shows a good agreement between our model results and his model results and his experimental results. For the 10” case some drift from the data can be observed. For longer spans (20”, 30” and 40”) our model tracks correctly with changes in span length. He also compared his model results with other test results. This was a relatively thin web, a 56 ga Polyester, with a Young’s Modulus of Ex = 658000 psi, an assumed Poisson’s ratio of ν = 0.3, and a width of W = 6”. The critical stress calculated from Eq.2.1 is 150 psi. The results from the Excel VBA code will be compared with both his model results and his experimental results below. 0 0.0005 0.001 0.0015 0.002 0 5 10 15 20 25 Mcr(in/in) Tension 56 ga Polyester, L=10" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.12 Wrinkles Due to Taper, 56 ga Polyester, L=10” 107 0.0002 0.0003 0.0008 0.0013 0.0018 0 10 20 30 Mcr(in/in) Tension 56 ga Polyester, L=20" Beisel's Model (COSMOS) Beisel Lab.Test Yurtcu's Model(VBA) Figure 5.13 Wrinkles Due to Taper, 56 ga Polyester, L=20” 0 0.0005 0.001 0.0015 0.002 0 5 10 15 20 25 Mcr(in/in) Tension 56 ga Polyester, L=30" Beisel's Model (COSMOS) B 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


