

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


CONNECTING WINDING CODES TO QUALITY CONTROL DEVICES By DIPESH DILIP MISTRY Bachelor of Science in Mechanical Engineering University of Pune Maharashtra, India 2007 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE July, 2010 iii CONNECTING WINDING CODES TO QUALITY CONTROL DEVICES Thesis Approved: Dr. J. K. Good Thesis Adviser Dr. R. P. Singh Dr. J. D. Jacob Dr. Mark E. Payton Dean of the Graduate College iv ACKNOWLEDGMENTS First of all, I would like to thank the Almighty for giving me the strength during the course of this research. I wish to express my sincere gratitude to Dr. J. K. Good, my adviser, for providing me an opportunity to work for him. He has been a great inspiration and motivated me all the way during my stay at OSU. I thank him for his wisdom and insight. He has been very considerate and supportive of me all the way. I would also like to thank my thesis committee members, Dr. R. P. Singh and Dr. J. D. Jacob. I wish to express my special appreciation to Mr. Ron Markum for guiding me through the experiments and helping me with anything I needed during this research. I thank my friends and colleagues who have always stood by me. Their opinions, suggestions and encouragement have been very helpful in various parts of my research. Finally, I am grateful to my parents, Mr. Dilip D. Mistry and Mrs. Darshana D. Mistry, and my family members, especially my elder brother, Mr. Bhupesh D. Mistry and my sister in law, Mrs. Upasana B. Mistry, who have provided immense psychological and moral support by constantly encouraging and believing in me throughout my life. Dipesh Dilip Mistry.v TABLE OF CONTENTS Chapter ........................................................................................................................... Page I. INTRODUCTION TO WINDING AND QUALITY CONTROL OF WOUND ROLLS ..... 1 II. LITERATURE REVIEW ................................................................................................ 6 2.1 Winding Models: .................................................................................................... 6 2.2 Wound Roll Contact Models: ................................................................................ 9 2.3 Quality Control Devices: .......................................................................................10 2.3.1 Schmidt Hammer: ......................................................................................... 11 2.3.2 PARO Tester: ................................................................................................ 11 2.3.3 Beloit Rho Meter: ..........................................................................................12 2.3.4 Online Density Measurement: ......................................................................12 2.4 Research Objective: ..............................................................................................13 2.5 Accomplishing the Objective: ...............................................................................14 III. EXPERIMENTAL PROCEDURE ................................................................................15 3.1 Winding Setup: .....................................................................................................15 3.2 Hardness Profile Measurement: .......................................................................... 17 3.3 Core Effect on Roll Hardness: ..............................................................................18 3.4 Hardness Measurements in MD and CMD: .........................................................19 3.5 Web Material Properties: .....................................................................................21vi 3.5.1 Stack Test: .....................................................................................................21 3.5.2 Stretch Test: ................................................................................................. 23 3.6 Properties of Polyester web used in this research (PET 377): ............................ 24 IV. DYNAMIC IMPACT MODEL ..................................................................................... 25 4.1 Model Development: ........................................................................................... 25 4.2 Finite Element Formulation: .............................................................................. 26 4.3 Modeling Rho Meter Parameters: ....................................................................... 32 4.4 Transient Analysis: ..............................................................................................40 4.5 Time Step Calculation: ........................................................................................ 42 V. RESULTS AND DISCUSSION .................................................................................... 45 5.1 Experimental Results: .......................................................................................... 45 5.2 Results from the Model: ...................................................................................... 53 5.3 Comparison of Results from Model with Experimental Results: ........................ 55 VI. CONCLUSION AND FUTURE SCOPE ......................................................................60 Future Scope: .............................................................................................................61 BIBLIOGRAPHY ............................................................................................................... 63 APPENDICES ................................................................................................................... 65vii LIST OF FIGURES Figure ............................................................................................................................. Page Figure 1.1: Quality Control Instruments for Wound Rolls ................................................. 4 Figure 3.1: High Speed Winding Line – Web Handling Research Center, OSU ...............16 Figure 3.2: Hardness Measurement of Wound Roll using the Rho Meter ........................ 17 Figure 3.3: Dog Bone Edge Formation due to Thickness Variation in the Web .............. 20 Figure 3.4: Experimental Setup for Stack Test ................................................................. 22 Figure 3.5: PET 377 Web Roll ........................................................................................... 24 Figure 4.1: Finite Element Formulation of Dynamic Impact Model ................................ 28 Figure 4.2: Internal Components of Rho Meter ............................................................... 33 Figure 4.3: Internal Components of Rho Meter (Top View) ............................................ 33 Figure 4.4: Drill Rods and Hinges Supporting the Moving Plastic Parts ......................... 34 Figure 4.5: Setup for Calculating Stiffness of Plastic Hinges ........................................... 35 Figure 4.6: Initial (A) and Final (B) Position of the Spring Mass System ....................... 36 Figure 4.7: Acceleration of the Rho Meter Striker Output on the Oscilloscope ............... 37 Figure 4.8: Initial Position of the Striker at Rest ............................................................. 39 Figure 4. 9: Finite Central Difference Estimates ..............................................................40 Figure 4. 10 Flowchart for Dynamic Impact Model .......................................................... 44 Figure 5.1: Roll Radius Vs Hardness at Constant Winding Tension ................................ 45viii Figure 5.2: Pressure Profile in Rolls Wound at Constant Tension ................................... 46 Figure 5.3: Comparison of Hardness of the Rolls Wound at Constant Tensions with Peak and Average Pressure in the Wound Roll ......................................................................... 47 Figure 5.4: Winding Tension Profile in the Roll due to Alternating Tensions ................. 49 Figure 5.5: Pressure Profile in Rolls Wound at Alternating Winding Tensions ............... 49 Figure 5.6: Comparison of Hardness Profile in the Roll Wound at Constant (5 PLI) and Alternating Tensions (5 PLI  3 PLI) ................................................................................ 50 Figure 5.7: Comparison of Hardness Profile in the Roll Wound at Constant (6 PLI) and Alternating Tensions (2 PLI  6 PLI) .................................................................................51 Figure 5.8: Hardness Comparison of Rolls in MD and CMD (PET 377) .......................... 52 Figure 5.9: Hardness Comparison of Rolls in MD and CMD (Mitsubishi Polyester) ...... 53 Figure 5. 10: Deceleration at Tip Node obtained by using different Time Steps (Δt) in Transient Analysis ............................................................................................................ 54 Figure 5.11: Comparison of Model Output for Linear and Non Linear Material Characteristics with Experimental data ........................................................................... 55 Figure 5.12: Dynamic Model Output for Different Mesh Size, Material Properties are Not Updated during Dynamic Analysis. .................................................................................. 56 Figure 5.13: Variation in Er calculated from Winding Pressures in an Element due to the Size of the Mesh. ................................................................................................................57 Figure 5.14: Dynamic Model Output for Different Mesh Size, Material Properties are Updated during Dynamic Analysis. .................................................................................. 59 ix LIST OF SYMBOLS AND ABBREVIATIONS WHRC Web Handling Research Center HSWL High Speed Web Line PET Polyester Film MD Machine Direction CMD Cross Machine Direction Er Radial Young‟s modulus Eθ Tangential Young‟s modulus εr Radial Strain εθ Tangential Strain σr Radial Stress σθ Tangential Stress σz Axial Stress νrθ Poisson‟s ratio of stress in θdirection to strain in rdirection r Radius of wound roll rc Outer radius of core x P Normal pressure K1 Material constant K2 Material constant C0, C1,..., Cn Constants in Hakiel‟s polynomial for pressure Tw Web tension t Web thickness ρweb Weight density of web material ρcore Weight density of core material h Web thickness made dimensionless by division by core radius δl Change in length te Thickness of quadrilateral element Elemental material matrix in Cartesian coordinates Elemental material matrix in cylindrical coordinates Transformation matrix θj Orientation for jth node [ke] Elemental stiffness matrix [me] Elemental mass matrix [KG], K Global stiffness matrix xi [MG], M Global mass matrix {F}, F (t) Global load vector [H] System matrix B StrainDisplacement matrix J Jacobian η, ξ Natural coordinates βk, βm Penalty number for stiffness and mass matrix respectively Keff Effective stiffness of drill rods Meff Effective mass of Rho meter striker and other components ωn First angular frequency of Rho meter striker d Displacement Velocity Acceleration Δt Time step for transient analysis Lmin Smallest element dimension in the mesh Cd Dilatational wave speed λ Lame‟s constant pli Pounds per linear inch 1 CHAPTER I INTRODUCTION TO WINDING AND QUALITY CONTROL OF WOUND ROLLS Webs can be defined as a material whose length to width ratio is very high and the thickness is very small as compared to both the length and width of the web. Webs are usually manufactured as thin, continuous and flexible layers of material. Webs compose a wide range of materials such as paper, plastic films, textiles, nonwoven materials, metal foils and composites. The most convenient form to store and transport webs is to wind them into rolls. These rolls are then unwound and rewound into subsequent production processes, which collectively adds value to these webs. Winding is the process by which the web is wound on a core to form a roll. The winding operation is done with machines called as winders. There are several types of winders used in the web handling industry. These are classified based on number of drums used. The most common types are single drum, multiple drum and beltreel winders. Based on the method by which the winding torque is applied, winders are categorized as center winders, surface winders, center winder with idler nips, center/surface winders, multiple nip winders, etc. 2 The process of winding is directly affected by the type of winder, the web material and the winder operating parameters (e.g. winding tension, speed, etc). The winding tension results in internal stresses that are responsible for giving the wound roll some integrity as a structure. The quality of the wound roll depends on the stresses which exist in it. Defects such as internal stress bursts, buckling, tearing, etc may occur in the wound roll if these stresses are not properly controlled. Thus, it is a stiff challenge for the web manufacturer and the winder operator to wind good quality rolls of their high quality products. The stresses in the wound roll are relatable to the roll defects. These stresses can be predicted by using winding models. Many winding models have been developed over the years for analyzing the stresses and pressures developed in the roll during winding. The majority of these winding models solve a second order differential equation in pressure written with respect to radius. This equation has non constant coefficients and must be solved numerically to produce results. Typical results consist of pressure and tangential stress profiles with respect to radius within the roll. Thus, winding models can be used to predict the optimum winding conditions for a given web material. Although, the winding model is a powerful quality control tool in the laboratory, it cannot determine the quality of the wound roll on the production floor. A winder operator must be able to check the wound rolls for defects on the shop floor. The most basic tool used for wound roll testing is “the backtenders friend,” a simple short length of stick used to manually strike the roll surface to judge the hardness across the roll width 3 by the sound and feel of the stick impacting the roll surface. With no means to quantify, the method limits the value of this test as a subjective determination of roll structure. (1) The most popular quality control devices on the shop floor are shown in Figure 1.1: Quality Control Instruments for Wound Rolls ,. These are handheld units such as the Beloit1 Rho Meter, the Schmidt2 Hammer and the PAROtester3. These devices measure the hardness of wound roll, which help in detecting the defects during winding. The Beloit Rho Meter, like the backtender‟s stick, is an impact tester. It measures the peak deceleration of a hammer striking the roll surface. The peak deceleration of the striker is converted to an arbitrary unit called “Rho”. The builtin scale permits quantifying the test results. The Schmidt Hammer (2) is another impact tester measuring the rebound of a plunger striking the roll surface. The magnitude of rebound is recorded on a builtin mechanical scale designated in “R” units. The PAROtest (3) is initiated by launching a spring loaded body against the test surface. The impact and rebound velocities are compared resulting in an instantaneous numerical hardness value designated in “L” units. This unit is related to the coefficient of restitution. 1 Beloit Manhattan Division. 1910 Lane Blvd. Kalamazoo, Michigan 49003. USA 2 Testing Machines, Inc. 2 Fleetwood Court, Ronkonkoma, NY 11779 USA 3 Testing Machines, Inc. 2 Fleetwood Court, Ronkonkoma, NY 11779 USA4 Figure 1.1: Quality Control Instruments for Wound Rolls (3), (2) The winding models output the stresses and pressures in the wound roll in terms of engineering units, such as pressure or stress. The quality control devices used on the shop floor output the roll hardness in arbitrary units of hardness. The arbitrary units of hardness cannot be directly coupled to the output in engineering units produced by the winding models. This study focuses on establishing a connection between the outputs of a winding model and a quality control device, the Rho Meter in this case. The Rho meter was selected due to availability in the labs of the Web Handling Research Center at OSU. Initially, experiments were conducted to study the nature of hardness variation in wound rolls as a function of the current outer radius of a winding roll. The relation between the peak pressure in the wound roll and the hardness of the wound roll was investigated. The dynamic hardness measurement range in terms of roll radius was also investigated. These experimental results were used to verify a model which combines a winding model and a dynamic impact model and can output the hardness of the wound roll in terms of 5 Rho units. This model can be an extension to all existing winding models. The experimental procedures and development of the model are discussed in subsequent chapters. The results obtained from the dynamic model are compared to the experimental results. 6 CHAPTER II LITERATURE REVIEW 2.1 Winding Models: Literature on the analysis of wound roll structure first appeared in 1960s. A wound roll structure can be defined as the web layers wound on a core at appropriate winding tension so that the wound roll has some integrity as a structure. Over the years, many winding models were developed for analyzing the stresses and pressures developed in the roll during winding. The early models appeared as 1D models with radial variation only. These models assumed a linear isotropic material and used stress formulas developed from the elasticity solutions of thick cylinders, in order to calculate the internal stresses of a wound roll. Although the 1D models could theoretically capture the accretive pressures which result from winding web layers onto a roll, they did not consider the material complexity due to anisotropy and radial nonlinearity that was required to produce the accurate results. Pfeiffer (4) investigated the radial nonlinearity of wound roll material. He found that the relation between the radial pressure and strain was exponential. 7 {2. 1} He developed a wound roll model in which the radial modulus was state dependent on pressure. {2. 2} Here K1 and K2 are material specific parameters calculated by fitting a curve between pressure versus strain data obtained from stack compression experiments of web material. Altmann (5) provided an analytical solution by assuming constant radial modulus, although different from circumferential modulus. Yagoda (6) was the first one to accurately treat the core boundary condition in wound rolls. Hakiel‟s work (7) is considered to be one of the most important steps in wound roll modeling. He used Altmann‟s formulation and developed a nonlinear anisotropic model. Hakiel incorporated his own model for the radial nonlinearity. He used a polynomial in pressure instead of Pfeiffer‟s exponential form. {2. 3} Hakiel‟s model was the first one to treat the state dependent properties and all boundary conditions rigorously. He derived a second order differential equation to calculate the increments in pressure ( ) in the wound roll. 8 {2. 4} He required two boundary conditions for the solution and so he utilized the outer and inner boundary conditions. The inner boundary condition is the compatibility of deformation of the core and the first layer of wound roll {2. 5} The outer boundary condition comes from the thin pressure vessel formula {2. 6} Hakiel employed a finite difference scheme to solve the set of equations for the addition of each lap and calculate the incremental pressure in each lap. The incremental pressures in each lap are summed with all previous increments to determine the total pressure. The material properties ( ) due to the current total pressure were then calculated and updated in the model prior to solving for the increments in pressure due to the addition of the next layer. After all layers are added the tangential stresses as a function of radius can be determined using the equilibrium equation posed in cylindrical coordinates: {2. 7} 9 2.2 Wound Roll Contact Models: Wound roll contact models have been developed over the years to analyze the contact problems that occur during winding or storage of wound roll. Ganapathi (8) developed a 2D axisymmetric finite element model to study the diametral compaction of wound roll, with state dependent properties, between rigid surfaces. He calculates the initial pressures in the wound roll by using Hakiel‟s (7) nonlinear winding model. The winding pressure is then used to calculate the radial modulus by using Pfeiffer‟s (4) relation. The wound roll cylinder is treated as an anisotropic cylinder for the formulation of the finite element contact model. The cylinder is meshed by using 4 noded quadrilateral elements (9). He adopted the general strategy of node to node contact in order to find the force which will bring down an associated node to the contacting surface. The contact is accomplished via linear interpolation wherein it requires three solutions per contacting node. The first solution is obtained by applying a unit load at the core and calculating the associated pull down force at the node. The next pull down force is obtained by applying half unit load and the final solution is obtained by calculating the load by linear interpolation of the first two results such that it will produce zero pull down force. At this point, the stresses are calculated and the material properties are updated. As a result of the contact pressures combining with the pressures due to winding, the radial modulus will now vary spatially with respect to radius and tangential location within the contact model. The operation is terminated when the total compression force is obtained. He verified his model experimentally (8). Mollamahmutoglu (10) developed a Nip and Wound Roll Contact model to analyze the contact pressures between nip and wound roll. The model is based on 2D plane stress finite element formulation of the contact of two orthotropic cylinders, one of 10 them represents the wound roll with nonlinear material properties and other represents the rigid nip roller. Similar to Ganapathi‟s model, the initial pressures due to winding are calculated by using Hakiel‟s (7) nonlinear winding model. The winding pressure is then used to calculate the radial modulus by using Pfeiffer‟s (4) relation. Mollamahmutoglu uses a very simple and efficient algorithm based on quasilinearization technique to obtain a rapid solution. The algorithm proposed by Mollamahmutoglu is much faster because it solves only once per contact node whereas the algorithm adopted by Ganapathi‟s model solves three times per node. This model was verified experimentally using Ganapathi‟s data and good comparison was found with Ganapathi‟s model when the problem of contact between wound roll and rigid wall was solved using the same technique. 2.3 Quality Control Devices: Smith (11) contends that roll hardness is probably the most important factor in determining the difference between a good and bad roll. Rolls that are wound too soft may have defects like starring and telescoping due to interlayer slippage during unwinding. Rolls that are wound too tight contain high inwound tension due to which tension bursts can occur inside the wound roll. The main aim of wound roll quality control measurement techniques is to determine if visually indiscernible defects are present into the roll. If the defects are detected then the rolls are rewound and future defects can be prevented by changing the parameters which affect the wound roll structure. Most of the quality control devices are handheld units and measure the roll hardness on the principle of impact. A few quality control 11 devices will be discussed briefly in this chapter. A thorough study of the Rho meter was conducted as we need to understand the instrument for modeling it into the dynamic impact finite element model that we intend to develop as a part of this research. This will be discussed in subsequent chapters. 2.3.1 Schmidt Hammer: The Schmidt test hammer is a handheld device based on the principle of rebound measurement. The Schmidt Hammer is composed of a spring loaded impact plunger which is pressed against the wound roll causing the compression of the spring. After reaching a certain point, the spring is released and an internal hammer mass is launched against the impact plunger. The amount of hammer mass rebound depends on the hardness of the wound roll. The magnitude of rebound is recorded on a builtin mechanical scale designated in “R” units. Wound roll hardness profile is generated by taking hardness measurements across the roll width. The hardness profiles are then used to detect the rolls with inconsistent hardness. In addition to the builtin scale, a paper tape print out is also available which gives an instant picture and records the cross roll hardness profile. (2), (1) 2.3.2 PARO Tester: PARO tester is an impact tester which is initiated by launching a spring loaded body against the wound roll surface. The impact and rebound velocities are measured and converted to output an instantaneous numerical hardness value. The value is the quotient of the impact and rebound velocities multiplied by 1000, results in output designated in “L” units. This unit is related to the coefficient of restitution. The accuracy and ease of implementation makes it one of the popular instruments on the shop floor. The digital display and inherent data memory makes it easy to interpret 12 and operate. The measurement information can be downloaded on a computer for further analysis. (3), (1) 2.3.3 Beloit Rho Meter: The Rho meter is an impact tester which measures the peak deceleration of a striker when it impacts on the roll surface. The striker is composed of steel and is a constant mass. A cantilever spring system with straight line motion accelerates the striker mass to a known constant velocity. A peak force is generated by the contact of striker with the surface of the wound roll. The peak force is detected by a high accuracy accelerometer, which is rigidly mounted on the striker to measure the peak deceleration on impact. The peak deceleration is converted into an arbitrary unit called “Rho”, which is a measure of roll hardness. The Rho meter has a digital display which indicates the roll hardness instantly and accurately. It is easy to operate and interpret. The Rho meter produces highly repeatable readings for a given roll. (4), (1) 2.3.4 Online Density Measurement: David Roisum‟s PhD work (12) was on online density measurement. The paper industry has been making online density measurements on winders for many years. It is easy to measure by using two encoders. One encoder is attached to an external nip and it measures the amount of web entering the wound roll and the peripheral velocity of the wound roll. A second encoder is attached to the core and measures the angular velocity of the core. During a measurement period one can determine 1. What length of web has entered the winder and, 2. How has the radius of the wound roll changed. 13 Given a web of known thickness one can determine its compressed volume and hence its density. Roisum tried to connect online density measurement to wound roll stresses through winding models. What he found at that time (1990) was that with encoder technology of that time he was marginally able to measure density with the accuracy he needed to predict wound roll deformations on materials as compressible as paper. Conceivably it would have worked better yet on materials with even greater compressibility like tissue. The stress condition inside the wound roll is a key factor in determining the quality of the roll. The winding process and associated parameters are generally responsible for the wound roll‟s stress condition. Stresses induced during winding are known to be related to defects. Therefore, determining the winding parameters that produce optimum roll stress profiles is of importance. Winding models can predict the stresses in terms of engineering units, which can be verified experimentally in a laboratory. On the production floor, the winder operator uses handheld quality control devices to detect the defects in the wound roll by measuring the roll hardness. 2.4 Research Objective: Review of literature shows that only one attempt has been previously made to provide outputs from wound roll models that can also be measured in the laboratory or on the production floor. That attempt was only marginally successful. The goal of this research is to develop extensions to wound roll models which converts outputs of pressures and stresses to units that can be measured using existing quality control devices, the Rho Meter in this case. 14 2.5 Accomplishing the Objective: The variation of roll hardness with respect to the roll geometry and winding conditions is investigated experimentally. The relation between the hardness and peak pressures in the wound roll is also studied. Experiments are also conducted to investigate the dynamic measurement range of the Rho Meter. The internal mechanism of the Rho Meter is studied and modeled into a finite element dynamic impact model with nonlinear material properties. The output from the finite element dynamic impact model is compared to the experimental results. The experimental procedure and development of dynamic impact model is discussed in detail in subsequent chapters. 15 CHAPTER III EXPERIMENTAL PROCEDURE In this research the hardness of wound rolls is measured using a Rho meter. Hardness variation of the rolls with respect to the outer radius is investigated. Rolls are wound at different tension levels. The winder is capable of holding the tension at zero velocity, thus facilitating the hardness measurements at intermediate outer roll radius positions from the core to the final wound roll radius. 3.1 Winding Setup: The High Speed Winding Line (HSWL) consists of number of rolls through which the web passes from unwind roller to the winding roller. It can reach high speeds of winding, by which the machine gets its name. The tension is applied using an automatic feedback tension control system. The winding machine is equipped with position guides which track and control the lateral position of the web during winding. The surface speed of the web is maintained low enough to avoid air entrainment problem. The webs were wound at five different levels of tension viz. 2, 3, 4, 5 and 6 PLI tension. 16 Figure 3.1: High Speed Winding Line – Web Handling Research Center, OSU 17 Figure 3.1 shows the HSWL. The rolls are wound on a steel core which has known properties and high stiffness. The procedure to measure the hardness of the roll is discussed in the following section. 3.2 Hardness Profile Measurement: Experiments were conducted on the web supplied by DuPont. The material was polyester 377 (PET 377) film in a 6inch width and 200 gage (0.002 inch) thickness. Since the HSWL is capable of holding tension at zero velocity, about half an inch (pile height) of web was wound on the steel core at a constant winding tension and then the machine was set to run at zero velocity. The Rho meter was used to measure the roll hardness at the current pile height of the web wound onto a core. Figure 3.2: Hardness Measurement of Wound Roll using the Rho Meter 18 Figure 3.2 shows the hardness measurement being taken by a Rho meter on a wound roll mounted on the HSWL. The hardness measurement was taken at five locations along the width of the web. The average hardness of the roll for the current pile height was calculated from these readings. The radius of the steel core was 1.7 inches and the final outer radius of the roll was 5.2 inches. The procedure was repeated for increments of one half inch of pile height. Thus, hardness readings were taken at seven radial locations in the wound roll. The hardness readings were plotted as a function of roll radius. The average hardness of the roll wound at different winding tension is compared to the peak pressure in the roll, which is calculated using a winding model. 3.3 Core Effect on Roll Hardness: The Rho meter makes dynamic measurement of hardness over a radial range. Experiments were conducted in order to investigate the range over which the hardness measurements are made by the Rho meter. The methodology adopted to study this range is by winding the web on a steel core (rigid) and then winding it on a comparatively soft core. The comparison of hardness at corresponding radial locations will give us an insight about the range of the Rho meter. We define “Core effect” as the variation in hardness of a roll due to winding it either on a steel (rigid) core or a softer core, at the same winding tension. The web was wound on a steel core at two different tension levels. These tensions were switched after about half inch increment in the pile height. The initial pile (0.5 inch), wound at low tension, acts as a soft core for the preceding pile (0.5 inch) of web that is wound at a higher tension. Experiments were conducted for two cases. For the 19 first case the higher and lower tension levels were 5 PLI and 3 PLI, respectively. For the second case, the tension levels were changed to 6 PLI for higher tension and 2 PLI for lower tension. In any case, as compared to the rigid steel core, each increment in pile height acts as a core for the preceding increment in pile to be wound onto the roll. Initially, the web was wound at high tension starting at the core and then the tension was reduced after one half inch pile height was wound on. The tension was alternated after every one half inch of pile height. After winding one half inch (pile height) of web, the HSWL is held at zero velocity while holding tension. The Rho meter was used to measure the roll hardness across the roll width and the average roll hardness was calculated at the current pile height of the web. Next, winding was started at low tension at the core and then the tension was increased after one half inch pile height. The tension was then alternated after every half inch of pile height. After winding one half inch (pile height) of web, the HSWL is held at zero velocity while holding tension. The Rho meter was used to measure the roll hardness across the roll width and the average roll hardness was calculated at the current pile height of the web. The final outer radius of the roll was 5.2 inches and hardness readings were taken at seven radial locations in the wound roll. The hardness readings were plotted as a function of roll radius. The hardness profile is compared to the pressure profile of the roll, generated in the winding model. 3.4 Hardness Measurements in MD and CMD: The striking mass of the Rho meter has a finite length. It is critical to determine the effect of contact length of the striker on the roll hardness measurements, when it 20 impacts the wound roll. The materials used for these experiments were the PET 377 previously described and a Mitsubishi polyester film in a 24inch width and 300 gage (0.003 inch) thickness. Mitsubishi intentionally produced this web with a highly radical thickness profile which results into very high variation in hardness across the width of the wound roll. Figure 3.3 shows the dog bone edge formation due to thickness variation in a Mitsubishi polyester wound roll. Figure 3.3: Dog Bone Edge Formation due to Thickness Variation in the Web The PET 377 web was wound on the HSWL at three different levels of tension viz. 1.5, 2 and 4 PLI. After winding about two inches (pile height) of web, the HSWL is held at zero velocity while holding tension. Now hardness measurements were taken by holding the Rho meter in the machine direction (MD) and the cross machine direction (CMD) at three locations along the width of the web. The Mitsubishi polyester web was 21 wound on HSWL at a tension of 1.5 PLI. After winding about ten inches (pile height) of web, the HSWL is held at zero velocity while holding tension. Again hardness measurements were taken by holding the Rho meter in MD and CMD at five locations along the width of the web. The hardness data obtained for both cases (MD and CMD) was plotted as a function of roll width. 3.5 Web Material Properties: 3.5.1 Stack Test: The stack test is conducted to measure the radial modulus (Er) of the web material. The radial modulus is a very critical variable in the winding process. Parameters like air entrainment, surface roughness and coatings combine to form the radial modulus of the wound roll, which is state dependent on the pressure in the wound roll. This nonconstant material property offers a challenge in understanding and modeling the wound roll structure. Figure 3.4 shows the experimental setup for stack test to measure the radial modulus (Er) of the web. 22 Figure 3.4: Experimental Setup for Stack Test The web is compressed slowly at constant velocity to allow any air to escape laterally form the stack as the compression progresses. Deformation and compression load are simultaneously recorded using a data acquisition system, during these tests. This data is used to infer the pressure as a function of strain in the stack. The experimental data is then used to produce a pressure versus strain chart. A least square approximation is used to fit a curve to the experimental data and the calculate K1 and K2 in Pfeiffer‟s exponential expression {2.1} for the pressure versus strain curve. K1 has the units of pressure and K2 is dimensionless. 23 3.5.2 Stretch Test: The stretch test is conducted to measure the tangential modulus (Eθ) of the web material. Tangential modulus of the web is an important factor in determining the pressures developed inside the wound roll. Both Er and Eθ affect the radial stresses and thus the pressure in the roll calculated using the second order differential equation {2.4}. The tangential modulus for any common engineering specimen could be measured using a material testing machine. For web materials better results are obtained by stretching a longitudinal length of the web (50 feet) and recording the load required to achieve a measured deformation. This method is preferred because the web must somehow be gripped to install a tension. Grip effects typically induce two dimensional stresses in the vicinity of the grips and thus the stress state in the web is not uniaxial. The long specimen minimizes the errors due to grip effects. In the stretch tests that are done 50 feet of web that is 6 inches wide are stretched. This results in an aspect ratio of 100, which serves to minimize the variability due to grip effects and specimen misalignment. One end of the web is fixed by taping it to the floor and a known load is applied at the other end. The load is constantly increased and the change in length ( δl ) is measured for application of each load. The loading is done for about twenty increments. The strain can then be calculated from the change in length for each load. The stress can be calculated from the area of the web and the force applied. The stress  strain curve is generated and the slope of this curve gives us the tangential modulus (Eθ). 24 3.6 Properties of Polyester web used in this research (PET 377): The 200 gage films have the following measured properties: Web width = 6 inches Web thickness = 0.002 inches Material constants for Pfeiffer‟s (4) relation calculated experimentally: K1 = 3.19 psi K2 = 38.35 Eθ = Ez = 711,000 psi Figure 3.5: PET 377 Web Roll 25 CHAPTER IV DYNAMIC IMPACT MODEL 4.1 Model Development: In this chapter a model is developed for a 2Dplane formulation of the dynamic impact of Rho meter striker and the wound roll. The model consists of two main components; a Winding model and a Dynamic Impact model. The winding model is a development similar to that of Hakiel (7) and the impact model is based upon a dynamic finite element formulation. Both models are developed using Excel and Visual Basic for Applications and the Excel/VBA code is shown in appendix A. The winding model employed in this analysis is based on Hakiel‟s non linear model for wound roll stresses (7). The input needed for this model are material properties which include the tangential modulus ( ) of the web, Young‟s modulus for the core, Pfeiffer‟s constants K1 and K2 (4) to establish the radial modulus of the web, Poisson‟s ratio for the web material and the core and the weight density of the web material and the core. The other input parameters include the geometry of the core and the roll, which include the inner and outer radii of the core, the outer radius of the roll, the thickness of the web and the width of the web material. The winding tension variation with radius is also required. 26 Winding formulations such as that proposed by Hakiel (7) solve for the increments in radial pressure within each layer of the roll due to the addition of the most recent layer. This is done by solving a second order differential equation written in terms of the incremental pressures which relies on two boundary conditions. One boundary condition is that the outside of the core and the inside of the roll must maintain continuity while the second condition assumes that the pressure beneath the outermost layer is known from the winding tension input and the current radius of the outer lap. Once the incremental pressures have been solved for the addition of the most recent layer they are summed with the incremental pressures within the laps due to the addition of all the previous layers to produce the total pressure within each lap. The total pressure in each lap is then used to update the radial direction modulus ( ) of the web using Pfeiffer‟s expression: , where P is the pressure due to winding. With updated properties a new layer is added and the second order differential equation is solved for the addition of the next layer until all layers have been added. After the final lap is added the final winding pressure and the radial modulus are known as functions of wound radius and will be used as initial values for the dynamic impact finite element model. 4.2 Finite Element Formulation: Figure 4.1 shows the finite element model for Rho meter striker impact with a wound roll. The finite element model is simplified by taking advantage of the symmetry and thus only a quarter of the roll is modeled. This simplification is justified because the striker impact effects are localized in the impact region only, which will be verified later, and due to the symmetry offered by the geometry of the wound roll. In this model each 27 node will be allowed to deform in the X and Y directions (u and v respectively), except for those nodes which may be constrained on the X and Y axes as shown in Figure 4.1. 2D quadrilateral plane stress elements with anisotropic properties are used in the formulation. Due to the fact that most of the deformation will occur near the impact region, a mesh refinement technique will be used along the radial and tangential direction that will densify the mesh in the impact region. The finite element formulation begins with mesh generation. The mesh generation consists of two parts. In the initial part a uniform mesh is generated in the radial direction for the core structure and in the second part a refined mesh is generated in the radial direction for the wound roll region. A refined mesh is generated in the tangential direction for the core and the wound roll region (10). 28 Figure 4.1: Finite Element Formulation of Dynamic Impact Model The initial wound roll material properties were computed by the winding model which determined Er as a function of r and E was assumed constant. These properties are now needed in the global Cartesian coordinates of the impact model. Thus the finite element formulation requires transformation of element material properties that are measured in the cylindrical coordinate (rθ) system of the winding model to the Cartesian coordinate (XY) system used in the finite element model. This transformation 29 is employed separately for each finite element during the calculation for the material matrix [ ]. {4. 1} Here is the orthogonal matrix which represents the transformation (10). For the element occupying a (i, j) position it can be given as: {4. 2} The stiffness matrix for each element [ke] is computed by using a 2x2 Gaussian quadrature scheme to numerically integrate (9): {4. 3} where J is the Jacobian and has the form: {4. 4} and, {4. 5} 30 {4. 6} {4. 7} {4. 8} These elemental stiffness matrices are assembled into a Global Stiffness Matrix [KG] by using the direct assembly procedure. The Global Stiffness Matrix [KG] is always symmetric and banded matrix. Although [KG] stored in a banded form saves some computer memory, in this case [KG] is stored as a square matrix in order to reduce the complexity of computations involved in the dynamic analysis. A lumped mass formulation was employed to compute the elemental mass matrix [me] for each element. The total element mass in each direction is distributed equally to the nodes of the element, and the masses are associated with translational degrees of freedom only (9). For the 2D quadrilateral element, the lumped element mass matrix [me] is: {4. 9} Further, the elemental mass matrices are assembled into a Global Lumped Mass Matrix [MG] by using the direct assembly procedure. The mesh refinement in the impact 31 region is very high and thus the lumped mass formulation is fairly accurate in this region. By using the lumped mass formulation a diagonal Global Lumped Mass Matrix [MG] is obtained, which offers computational efficiency in the dynamic analysis algorithm that will be discussed later. Geometric constraints must be enforced to the finite element model, which occur due to symmetry. As seen in figure 4.1, we need to constraint the vertical movement of the nodes on the lower surface (Y=0) and the horizontal movement on the left boundary (X=0) of the quarter cylinder. Initially, the penalty method was used to enforce these geometric boundary conditions. A penalty number (βk) is added to the appropriate components of the global stiffness matrix [KG] and an inertial penalty number (βm) is added to the appropriate components of the global mass matrix [MG]. The penalty approach resulted in huge numbers on the major diagonal of the global stiffness matrix, which required the value of critical time step (Δt) for the transient analysis to be reduced accordingly. Thus the overall time period over which the deceleration of the striker mass is calculated, reduced significantly and the solution never reached the peak deceleration. Due to the difficulties encountered by the penalty approach, the elimination approach was employed for enforcing the geometric constraints. In this method, the constraints are enforced by eliminating the rows and columns in the global stiffness matrix [KG] and the global mass matrix [MG], corresponding to the constrained degree of freedom. After enforcing the constraints, we obtain a reduced global stiffness and mass matrix. The advantage of the elimination approach over the penalty approach is that the value of critical time step (Δt) is increased resulting in a longer time step for transient analysis and decreased computational time. The Rho meter dynamic elements must also be included into the finite element model which is discussed in the following section. 32 4.3 Modeling Rho Meter Parameters: The internal mechanical components of a Rho meter are shown in figures 4.2 and 4.3 below. The striker impacts the wound roll at a known velocity and accelerometer measures the peak deceleration of the striker after the impact. The striker has a straight line motion which is accelerated by a cantilevered spring system to a desired contact velocity. The most important parameters that affect the peak deceleration on impact are the effective stiffness of the cantilevered spring system, the effective mass of the system, the initial velocity of striker and the hardness of the wound roll. For simplicity, the entire Rho meter is not modeled. The important Rho meter parameters are included at the tip node location as shown in Figure 4.1. Due to this, there are no external forces acting on the system and thus the global force vector {FG} is always equal to zero through time in the dynamic model. 33 Figure 4.2: Internal Components of Rho Meter Figure 4.3: Internal Components of Rho Meter (Top View) 34 The effective stiffness of the cantilevered spring system comes from the stiffness of the two steel drill rods and the stiffness of the composite hinges inserted in the moving plastic parts. These components are shown in Figure 4.4. The two drill rods measure 4.75" in length and 0.095" in diameter. These rods are made from steel. With known dimensions and material properties, the stiffness of each of the two drill rods was calculated from equation 4.10; {4. 10} Figure 4.4: Drill Rods and Hinges Supporting the Moving Plastic Parts 35 The stiffness of the hinges supporting the moving plastic parts was determined experimentally by measuring the displacement of the spring mass system and the corresponding force by using a handheld force gage. Figure 4.5 shows the setup for measuring the force and the corresponding displacement of the spring mass system. Figure 4.5: Setup for Calculating Stiffness of Plastic Hinges 36 The stiffness of the drill rods come into consideration after the spring mass system is displaced by 0.0625 inches. This is a nonstandard configuration for the Rho meter due to some rubber grommets that were missing in one of the WHRC‟s two units. The force applied by a handheld force gage up to a displacement of 0.0625 inches is directly affected by the stiffness of the hinges supporting the moving plastic mass. Figure 4.6 shows the initial and final positions of the spring mass system before and after the force is applied by the force gage. Figure 4.6: Initial (A) and Final (B) Position of the Spring Mass System The force required for producing a deformation of 0.0625 inches was found to be 0.58 lb. The stiffness of the hinges was calculated by using these readings; {4. 11} Thus, the effective stiffness of the system will be the sum of the stiffnesses of the drill rods and the plastic hinges. 37 {4. 12} The striker of the Rho meter is expected to stay in contact with the wound roll surface during the deceleration period. Thus, the effective stiffness ( ) of the system is simply added to the tip mode location of the global stiffness matrix [KG]. The natural angular frequency ( ) of the system can be measured experimentally by connecting the output of the accelerometer to an oscilloscope. The acceleration sensed by the accelerometer is output on the oscilloscope to measure the natural frequency of the system. The natural frequency of the system is calculated from the two consecutive peak acceleration observed after the impact of the striker. Figure 4.7 shows the acceleration of the Rho meter striker output on the oscilloscope. Figure 4.7: Acceleration of the Rho Meter Striker Output on the Oscilloscope 38 The natural angular frequency ( ) of the system observed on the oscilloscope was 23.80 Hz. Thus; {4. 13} The effective mass ( ) of the system is a summation of the mass of the striker, the support rod, the aluminum spring support, the moving plastic mass, the accelerometer, and the effective mass of the cantilevered drill rods shown in figures 4.2 and 4.3. The effective mass ( ) of the system can be calculated from equation; {4. 14} The effective mass ( ) of the system is simply added to the tip mode location of the global lumped mass matrix [MG]. The initial velocity of the Rho meter striker can be calculated from the first natural angular frequency ( ) of the system and assuming a harmonic displacement function. A harmonic displacement function (d) was assumed of the form: {4. 15} {4. 16} when the striker is released, ; where do is the height at which the striker is released above the wound roll by the Rho meter instrument. When is , the velocity will be maximum and will be zero. 39 {4. 17} Figure 4.8 shows the initial position of the striker at rest. As observed in Figure 4.8, the striker is not in line with the surface of the Rho meter shoe plate at time (t) =0. The distance by which the striker needs to displace in order to strike the surface of the wound roll was measured to be 0.125 inches. The time required for this displacement of the striker was calculated to be 7.016 x 103 sec. Figure 4.8: Initial Position of the Striker at Rest Thus the velocity of the striker on impact with the surface of the wound roll will be; 40 {4. 18} The velocity of the striker on impact ( ) is added to the tip node location of the global velocity vector at the start of the transient analysis. With the Rho meter parameters calculated and input within the appropriate locations in the global stiffness matrix, global mass matrix and global velocity vector, the transient analysis can be performed which is discussed in the following section. 4.4 Transient Analysis: The Direct Integration method is used to perform the transient analysis of the impact of Rho meter striker and the wound roll. In this method the derivatives of the deformation through time are evaluated by using Finite Central Difference method. The central difference estimates of the velocity and acceleration vectors at time ti are calculated using the following relations; Figure 4. 9: Finite Central Difference Estimates 41 {4. 19} {4. 20} The calculation of the time step for the transient analysis is very critical because it directly affects the accuracy of the solution. The calculation of time step will be discussed in the subsequent section. In the dynamic code, the peak deceleration is sought after the impact. The procedure applied in the dynamic code for calculating the decelerations is as follows: 1. Calculate time step 2. Read initial conditions: , which are known. 3. Solve for previous time step; {4. 21} 4. With determined, now solve for ; (for next ) {4. 22} This equation results from inserting 4.20 into the equation of motion shown in 4.22, (13): {4. 23} 5. With given and determined, we now solve for {4. 24} 6. Calculate corresponding deceleration using {4.15} or : 42 {4. 25} 7. Use steps 56 in a recursive fashion until the first peak deceleration is obtained. Once the peak deceleration of the striker is detected by using the above procedure, it is converted into hardness in terms of “Rho‟s”. The conversion is obtained by using the following relation that was used by the manufacturer of the Rho meter: {4. 26} The result, in terms of Rho hardness, is then printed on the Excel sheet. 4.5 Time Step Calculation: The explicit dynamic finite element analysis is very sensitive to the selection of the time step . The explicit procedure integrates through time by using small time increments as discussed in the preceding section. The value of Δt is very critical, because it affects the results directly. If too high a value for Δt is entered, divergent output results. If too low a value for Δt is entered, large rounding errors can result. The time step is generally expressed in terms of the highest natural frequency of the system as given in the following equation; {4. 27} The highest natural frequency ( ) can be calculated by using vector iteration method. In this method, if we use the system matrix , then the solution is converges to the highest natural frequency. The iteration process requires huge 43 computational times for the calculation of time step. Thus we switched to a more efficient method discussed below. A second alternative method allows the time step to be approximated at the elemental level. The time step ( ) for the dynamic problem is calculated based upon the smallest transit time of a dilatational wave across any of the elements in the mesh (14). This results in saving the computational time associated with the vector iteration method. The approximate value for can be calculated from equation {4.28}; {4. 28} Where, {4. 29} „Lmin‟ is the smallest element dimension in the mesh and „Cd‟ is the dilatational wave speed expressed in terms of Lame‟s constant and material density. This estimate for Δt is only approximate and is typically scaled by a factor between 1 and for two dimensional analysis and between 1 and for three dimensional analysis (14). For this case the value obtained by the method discussed above needs to be further scaled by a factor between 0 and 1 in order to compensate for the material non linearity parameters. The value of this factor is dependent on the mesh density of the wound roll. A flow chart of the dynamic impact finite element model is shown in figure 4.10.44 If Dtip=1st Peak Dec YES A Apply geometrical constraints Begin Transient Analysis. 1) Calculate Δt 2) Read Initial Velocity 3) Form initial deformation vector Use Finite – Central Difference Method to estimate the deformations. Calculate acceleration from deformations at previous time steps Calculate deformations at tip node for next Δt Calculate corresponding peak deceleration (Dtip) at tip node Calculate elemental stresses and corresponding radial pressures. Update Global Stiffness matrix. Apply Rho meter parameters & Constraints Calculate Rho Hardness Print Output =Rho Hardness NO STOP START Input: Material Properties, winding parameters, core and roll geometry. Execute Hakiel type Winding model and calculate Pr & Er for the wound roll. Generate Finite element mesh, using mesh refinement in the web region Generate elemental stiffness matrices using Er due to winding. Assemble Global stiffness matrix. Generate elemental mass matrices. Assemble Global Lumped mass matrix. Model Rho meter parameters: Add effective stiffness of the cantilevered spring system to tip node location in the Global Stiffness matrix. Add effective mass of striker to tip node location in the Global Mass matrix. A Figure 4. 10 Flowchart for Dynamic Impact Model 45 CHAPTER V RESULTS AND DISCUSSION 5.1 Experimental Results: The PET 377 web was wound at different winding tensions viz. 2 PLI to 6 PLI. The Rho meter was used to measure the hardness of the roll across the width. The average roll hardness was calculated at the current pile height (radius) of the web. The measurements were taken at seven radial locations. Figure 5.1 shows the hardness profile for the rolls wound at different tensions. Figure 5.1: Roll Radius Vs Hardness at Constant Winding Tension 0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 6 Hardness (Rho) Roll Radius (inch) 2 PLI 3 PLI 4 PLI 5 PLI 6 PLI Winding Tension46 The hardness values obtained at lower winding tensions is low and it increases as the winding tension is increased. The hardness near the core is higher as compared to the hardness away from the core, although, this trend becomes less pronounced at higher winding tensions. It was observed that at low winding tensions (23 PLI), the difference in hardness between consecutive winding tensions is greater as compared to the difference in hardness of rolls wound at higher tensions (46 PLI). Thus it can be inferred that there is a nonlinear relation between the wound roll hardness and the winding tensions. The hardness profile was compared to the pressure profile of the wound roll. Figure 5.2 shows the pressure profile for the rolls wound at different tensions. Figure 5.2: Pressure Profile in Rolls Wound at Constant Tension 0 50 100 150 200 250 300 350 400 450 500 0 1 2 3 4 5 6 Presssure (psi) Roll Radius (Inch) 2 PLI 3 PLI 4 PLI 5 PLI 6 PLI Winding Tension47 The pressures inside the wound roll were generated by using Hakiel‟s winding model. The pressure profile obtained from the winding model is of a plateaued roll where the pressures may be nearly constant throughout the roll. It was observed that, unlike the hardness, the pressure increases exponentially with an increase in winding tension. Figure 5.3: Comparison of Hardness of the Rolls Wound at Constant Tensions with Peak and Average Pressure in the Wound Roll Figure 5.3 shows comparison of the hardness of rolls wound at constant tension and the peak and average pressures in the roll. The pressure in the roll is generated by using Hakiel‟s winding model. The peak pressure occurs at the core of the wound roll, whereas the average pressure is calculated from the pressures at different radial locations within the wound roll. It was observed that the peak pressure increases 0 10 20 30 40 50 60 70 80 0 50 100 150 200 250 300 350 400 450 500 0 1 2 3 4 5 6 7 Hardness (Rho) Pressure (psi) Winding Tension (PLI) Peak Pressure Average Pressure Hardness: Experiments Roll Outer Radius = 5.2 inch48 exponentially with an increase in winding tension. The average pressure in a wound roll is comparatively lower than the peak pressure, but shows a similar trend with an increase in winding tension. The roll hardness shows a different trend. The slope of the hardness curve decreases with an increase in winding tension. Thus it can inferred that after a certain tension level, the hardness measurements obtained by the Rho meter will reach an asymptote. As observed in Figure 5.1, the hardness readings of the roll near the core are higher due to the effect of the core. The core in this case is radially much stiffer than the web layers in compression. This core effect is observed because the Rho meter makes dynamic measurement of the roll hardness over a localized range. Rolls were wound by alternating the winding tension after every one half inch pile height and the hardness measurements are taken across the roll width. Figure 5.4 shows the winding tension profile in the rolls wound at alternating tension. The pressure profile for these rolls was generated using Hakiel‟s winding model and is shown Figure 5.5. Experiments were conducted for two cases. For the first case the higher and lower tension levels were 5 PLI and 3 PLI, respectively. For the second case, the tension levels were changed to 6 PLI for higher tension and 2 PLI for lower tension. 49 Figure 5.4: Winding Tension Profile in the Roll due to Alternating Tensions Figure 5.5: Pressure Profile in Rolls Wound at Alternating Winding Tensions 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 Winding Tension (PLI) Roll Radius (inch) Alternating Tensions 5 PLI 3 PLI Alternating Tensions 2 PLI 6 PLI STEEL CORE Web Roll 0 30 60 90 120 150 180 210 240 270 300 330 360 0 1 2 3 4 5 6 Pressure (psi) Roll Radius (Inch) Alternating Tension 5 PLI 3 PLI Alternating Tension 2 PLI 6 PLI50 The hardness profiles for the rolls wound at alternating tension are shown in Figure 5.6 and Figure 5.7. It was observed that the average hardness of the roll wound on a steel core is higher when compared to the average hardness of the rolls wound on a core, whose stiffness is lower than the steel core. Thus, it can be inferred that the Rho meter hardness measurements are sensitive to changes in winding tension. With the data available from this experiment, it is difficult to find the radial range over which the Rho meter makes the hardness measurements. In order to investigate the “Core Effect” and the radial range of the Rho meter, the hardness measurements need to be made at relatively lower pile heights than half of an inch. These measurements would be effective in capturing the actual trend of the roll hardness profile in the radial direction. The experimental results can be used to study the variation in hardness of the wound roll due to stiffness of the core. Figure 5.6: Comparison of Hardness Profile in the Roll Wound at Constant (5 PLI) and Alternating Tensions (5 PLI  3 PLI) 0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 6 Rho Hardness (Rho) Roll Radius (inch) Alternating Tensions 5 PLI 3 PLI 5 PLI Constant Tension51 Figure 5.7: Comparison of Hardness Profile in the Roll Wound at Constant (6 PLI) and Alternating Tensions (2 PLI  6 PLI) Experiments were conducted to determine the effect of contact length of the Rho meter striker on the roll hardness measurements, when it impacts the wound roll. Two different web materials were used for this experiment, PET 377 and Mitsubishi polyester. PET 377 web was wound at different winding tensions and hardness readings were taken across the width of the web for each winding tension. The hardness measurements were taken by holding the Rho meter in the machine direction (MD) and the cross machine direction (CMD) at these locations. 0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 6 Rho Hardness (Rho) Roll Radius (inch) Alternating Tensions 2 PLI 6 PLI 6 PLI Constant Tension52 Figure 5.8: Hardness Comparison of Rolls in MD and CMD (PET 377) Figure 5.8 and 5.9 show the comparison of hardness of the roll for measurements taken in MD and CMD at the same locations along the width. In Figure 5.9 it can be observed that the hardness of the wound roll is very high at one of the edges. This results from the formation of dog bone edge due to very high thickness variation in the web. 0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 Hardness (Rho) Roll width (inch) 1.5 PLI MD 2 PLI MD 4 PLI MD 1.5 PLI CMD 4 PLI CMD 2 PLI CMD Winding Tension Roll OuterRadius = 4.2 inches53 Figure 5.9: Hardness Comparison of Rolls in MD and CMD (Mitsubishi Polyester) Good comparison was observed for the hardness measurements taken in MD and CMD. Thus it can be inferred that the contact length of the striker on impact has a less pronounced effect on the hardness measurements. 5.2 Results from the Model: The winding model and the finite element dynamic model is developed using Excel and Visual Basic for Applications. Initially, the convergence of the model was verified by comparing the outputs obtained for a constant set of input data by varying the time step (Δt) for transient analysis. Three different values for Δt were used for the transient analysis and the deceleration at the tip node was plotted as a function of time 0 20 40 60 80 100 120 140 160 0 5 10 15 20 25 Hardness (Rho) Roll width (inch) 1.5 PLI MD 1.5 PLI CMD Winding Tension Roll OuterRadius = 10 inches54 for all the three cases. Figure 5. 10 shows the deceleration at the tip node obtained by using three different values for Δt in the transient analysis. The deceleration curve from the model has a similar trend and the peak deceleration occurs at the same time for all the three cases. Thus it can concluded that the model shows good convergence when the given time step (Δt) is less than the critical time step for the problem. Figure 5. 10: Deceleration at Tip Node obtained by using different Time Steps (Δt) in Transient Analysis 55 5.3 Comparison of Results from Model with Experimental Results: The model output was compared to the experimental data. The non linear material properties were employed into the model and the material properties were updated after each solution in time step. Figure 5.11: Comparison of Model Output for Linear and Non Linear Material Characteristics with Experimental data Figure 5.11 shows the comparison of model output with the experimental data for two different rolls wound at constant tensions of 2Pli and 6 Pli respectively. The hardness output in terms of Rho units is obtained from the model at different outer roll radius locations. The output obtained from the model for lower winding tension (2 Pli) is 0 10 20 30 40 50 60 70 0 1 2 3 4 5 6 Hardness (Rho) Roll Radius (inch) Experiment 2 PLI Non Linear (Model) 2 PLI Experiment 6 PLI Non Linear(Model) 6 PLI56 in better comparison with the experimental data as compared to the model output for higher winding tension (6 Pli). This difference may be due to the coarse mesh size used for the above simulations. During the dynamic analysis the material properties are updated after each solution in time step, which results in huge computational times. The coarse mesh size was used in order to reduce the computational time for dynamic analysis. Figure 5.12: Dynamic Model Output for Different Mesh Size, Material Properties are Not Updated during Dynamic Analysis. 57 Figure 5.12 shows the model output for different mesh sizes. In order to reduce the computation times, the material properties are not updated during the dynamic analysis. As the mesh size increases, the deceleration output from the model decreases. Typically in FEA increased mesh size leads to results that are more accurate but in fact this is not true in this case. This is due to the estimation of the radial modulus (Er) of the elements in the wound roll region is calculated from the average pressures obtained from the winding model. Figure 5.13: Variation in Er calculated from Winding Pressures in an Element due to the Size of the Mesh. Figure 5.13 shows the variation in Er, calculated from the winding pressures, due to the size of the element in the mesh. As the mesh size gets coarser, the size of the element in the mesh gets larger and thus the average pressure within the element increases resulting in a higher value of Er for the element. Thus, the stiffness of the elements in the coarse mesh is higher as compared to the stiffness of the elements in the fine mesh. The 58 higher stiffness of the elements in the web region results in higher values for deceleration of the Rho meter striker during the dynamic analysis. Thus, if the material properties are not updated during the dynamic analysis, the output obtained from the model for different mesh sizes will not be accurate and huge errors result due to neglecting the change in material properties due to localized increase in pressures during the impact between the Rho meter striker and the wound roll. Next, the output from the dynamic model was obtained for different mesh sizes by updating the material properties after each solution in time step during the dynamic analysis. Due to this the computation times were high for this model and thus the deceleration output from the model was obtained only for a few time steps. Figure 5.14 shows the model output for different mesh sizes. As the mesh size increases, the deceleration output from the model increases. This will produce model results in Figure 5.11 that are closer to the test results. This is due to the fact that the Er for the elements in the web region is a function of the total pressure (PrTotal) due to winding and the localized pressures due to impact of the Rho meter striker. 59 Figure 5.14: Dynamic Model Output for Different Mesh Size, Material Properties are Updated during Dynamic Analysis. As the element size becomes smaller, the average pressures calculated at the four Gauss points within the element will be more accurate. The total pressure (PrTotal) in the element will result into higher Er for the element and thus the stiffness of the element will increase with a finer mesh size. These results support the argument that a finer mesh sizes will produce more accurate results from the model. 60 CHAPTER VI CONCLUSION AND FUTURE SCOPE In this research the experimental data was utilized to develop a model for dynamic impact of the Rho meter striker with the surface of the wound roll. From the results obtained from the experiments we can conclude that; The hardness variation in the wound roll with respect to the winding tension and the outer roll radius is studied experimentally. The hardness of the wound roll increases with an increase in the winding tension. Pressures and Hardness increase with an increase in tension, but the slope of the curves show a different trend. It can be inferred that the hardness measurements obtained by the Rho meter will increase up to a certain point, with an increase in winding tension, after that the measurements obtained by the Rho meter will reach an asymptote. The Rho meter makes a dynamic measurement of roll hardness over a range of wound roll radius close to the outer surface where the impact occurred. The output from the device is sensitive to variations in winding tension and the stiffness of the core. Further experimentation is required to study what portion of the outside of the winding roll affect the measurements that are made by the Rho meter. 61 The hardness measurements made by the Rho meter by aligning the axis of the Rho meter striker in MD and CMD are similar. Thus the geometry of impact of the Rho meter striker against the wound roll does not affect the hardness measurements by a considerable amount. The geometric constraints enforced by using the Elimination approach works better than the geometric constraints enforced by using the Penalty approach for this case. Although, the elimination approach increases the complexity due to change in dimension of the related matrices, they result into lower computation times. The algorithm used in the model for calculating the time step (Δt) shows good convergence for a given set of inputs. The selection of Δt is influenced greatly by the mesh density. The assumptions made in the development of the dynamic impact model are reasonable as the model results are in fair agreement with the experimental results. The output obtained from the dynamic model is in better comparison with the experimental data only if the material properties are updated locally during the dynamic analysis. Although the time required for computations is very high in this case the properties must be updated to produce accurate results. The output from the model obtained by employing a fine mesh size is more accurate as compared to the output obtained by employing a coarse mesh size. Future Scope: The main idea behind this study was to develop an extension to the winding models which can convert the outputs of pressures and stresses to units that can be measured using existing quality control devices. We have successfully modeled the Rho 62 meter in this research but it has been shown that reasonable mesh densities and thus large computation times are required to produce good converged results. To this end the computational times of the model can be reduced by using a localized impact model instead of using a quarter model of the wound roll, as modeled in this case. The region affected by the impact can be determined by examining the variation in pressures of the elements. The region where the pressure variation due to impact is negligible can be omitted and only the affected region could be modeled. This will result in more accurate results with reduced computational times for the dynamic analysis. To get more confidence in this model, experiments similar to the ones we conducted in this study need to be conducted for different web material and the results should be compared with the output from the model. This work could be extended further to develop similar models for some other popular quality control devices like the Schmidt hammer and the PARO tester. 63 BIBLIOGRAPHY 1. Wound Roll Testers. Paper Industry Web (PIW). [Online] www.paperindustryweb.com/rolltest.htm. 2. Paper Roll Hardness Tester. Schmidt Roll Hardness Tester. [Online] Testing Machines Inc. http://www.testingmachines.com/4020paperrollhardnesstesterrecording.html. 3. Paper Roll Hardness Tester. PARO Tester2. [Online] Testing Machine Inc. http://www.testingmachines.com/4022parotester2.html. 4. Pfeiffer, J. D., Internal Pressures in a Wound Roll of Paper. 5, August 1966, TAPPI Journal, Vol. 49, pp. 342  347. 5. Altmann, H. C., Formulas for Computing the Stresses in CenterWound Rolls. 4, April 1968, TAPPI Journal, Vol. 51, pp. 176  179. 6. Yagoda, H. P., Resolution of Core Problem in Wound Rolls. 4, Transactions of ASME, December 1980, Journal of Applied Mechanics, Vol. 47, pp. 847  854. 7. Hakiel, Z., Nonlinear Model for Wound Roll Stresses. 5, May 1987, TAPPI Journal, vol 70, no. 5, Vol. 70, pp. 113  117. 8. Ganapathi, S., Diametral Compression of Wound Rolls with State Dependent Properties. MS Thesis, Oklahoma State University, December 2007.64 9. Chandrupatla, T. R. and Belegundu, A. D., Introduction to Finite Elements in Engineering. 3. Prentice Hall, Inc, 2003. 10. Mollamahmutoglu, C. A 2D AxisSymmetric Wound Roll Model Including Nip Effects. PhD Dissertation, Oklahoma State University, December 2009. 11. Smith, R. D., The Art of Winding Quality Rolls. August 2001, Paper Film & Foil Converter, pp. 46  53. 12. Roisum, D. R., The Measurement of Roll Stresses During Winding. PhD Dissertation, Oklahoma State University, 1990. 13. Komzsik, L., Computational Techniques of Finite Element Analysis. CRC Press, 2005. pp. 195  202. 14. Kandadai, B., The Development of WoundOnTension in Webs Wound into Rolls. PhD Dissertation, Oklahoma State University, December 2006. APPENDICES APPENDIX A The dynamic impact model developed using Excel and Visual Basic for Applications is shown below: '**********************************************************************************'‟*********DYNAMIC IMPACT MODEL FOR CONNECTING WINDING CODES******** ‟*************************TO QUALITY CONTROL DEVICES************************ '**********************************************************************************' Option Explicit '**************INPUT PARAMETERS Dim nw, nc, nt As Integer Dim rcorein, rwebin, rwebout, te As Double Dim Erweb, Etweb, vrtweb, Ercore, Etcore, vrtcore, rhoweb, rhocore As Double Dim K1, K2 As Double Dim alpha, beta1, MultiFactor As Double Dim start '**************VARIABLES FOR MESH REFINEMENT AND OTHER CALCULATIONS Dim stepfactor, powerfactor, traystep, hstep, QWEBC(), RWEBC() As Double Dim DTHETA, hweb, hcore, Pi, chmax, chmin As Double Dim PR(), PRI(), um(), PRTotal(), PRStress() As Double Dim XC(), YC(), radpos As Double Dim THETAC() As Double Dim p(8) As Integer Dim gp(2) As Double Dim NN As Integer Dim KG(), FG(), Kdummy(), Kconstrained(), Mdummy(), Mconstrained() As Double Dim mpa, x1, x2, x3, x4, y1, y2, y3, y4, RD(3, 3) As Double Dim xi, xo, yi, yo, ri, ro, THETAE As Double Dim J11, J12, J21, J22, BL(3, 8), matcon As Double Dim ke(8, 8), D(3, 3) As DoubleDim t, s, detJ As Double Dim sf1t, sf2t, sf3t, sf4t, sf1s, sf2s, sf3s, sf4s As Double Dim q1, q2, q3, m1, m2, m3, m4, z1, z2, n1, t1, t2, t3 As Integer Dim MG(), CG(), MGINV() As Double Dim Mele(8, 8), NTN(8, 8) As Double Dim PC As Double Dim mpc, meff, keff, Xconstraints, Yconstraints, XCons, YCons As Integer Dim Irod, Drod, Erod, Lrod, Krod As Double '**************WINDER VARIABLES Dim KW(), AW(), BW(), CW(), PW(), PIW(), tw(), RW(), ERW() As Double Dim beta, FR, TWI, TWF, CS, hWINDER As Double Dim NL As Integer '**************DYNAMIC PROBLEM VARIABLES Dim delT, cons, Rhos As Double Dim U00(), U0(), U1(), Un(), Udummy(), Ucopy(), H() As Double Dim V00(), V0(), V1(), Vn() As Double Dim A00(), A0(), A1(), An(), Ag() As Double Dim F0(), MinF(), MinM(), KUC(), G() As Double Dim invel, tipnode, vec1, nmax As Integer Dim NE As Integer Dim cstress(3) As Double Dim HU0(), HUn(), MinMU0(), MinMU00() As Double Dim C, Omega, Iter1(), Iter2(), delTmax, delTnl As Double Dim Lmin, Cd, lamda, mu As Double Dim MTI2(), HmDelT(), BrTI(), BrTDI(), BrMinv(), M2d0(), KT2d0(), KT2() As Double '**************CHARTS 'Dim LoaddefChart, PressureChart, WPressureChart, NCZChart, MPCChart, NodeChart As Object Dim length As Double '********************************************************************************** '********************************************************************************** Sub Main() start = Timer Call Inputdata Application.ScreenUpdating = False Application.DisplayAlerts = False Call CartesianMesh Call RollWinder Call StiffnessSystem Call MassSystemLumped Call MassEffective Call StiffnessConstraint Call MassConstraint Call MassLumpedInverse Call TimestepAbaqus Call TransientDI 'Call CartesianResults Worksheets("INPUT").Activate Range(Cells(31, 2), Cells(31, 2)) = Timer  start End Sub '********************************************************************************** Sub Inputdata() 'mesh parameters and geometrical data hWINDER = Range(Cells(13, 2), Cells(13, 2)) nw = Range(Cells(4, 2), Cells(4, 2)) nt = Range(Cells(5, 2), Cells(5, 2)) nc = Range(Cells(6, 2), Cells(6, 2)) rcorein = Range(Cells(8, 2), Cells(8, 2)) rwebin = Range(Cells(9, 2), Cells(9, 2)) rwebout = Range(Cells(10, 2), Cells(10, 2)) te = Range(Cells(11, 2), Cells(11, 2)) 'material properties Etweb = Range(Cells(16, 2), Cells(16, 2)) vrtweb = Range(Cells(17, 2), Cells(17, 2)) Ercore = Range(Cells(19, 2), Cells(19, 2)) Etcore = Range(Cells(20, 2), Cells(20, 2)) vrtcore = Range(Cells(21, 2), Cells(21, 2)) rhocore = Range(Cells(26, 2), Cells(26, 2)) rhoweb = Range(Cells(27, 2), Cells(27, 2)) MultiFactor = Range(Cells(28, 2), Cells(28, 2)) K1 = Range(Cells(14, 2), Cells(14, 2)) K2 = Range(Cells(15, 2), Cells(15, 2)) 'Pi constant Pi = 3.141592654 'Calculate layer thickness for core and web regions hweb = (rwebout  rwebin) / nw hcore = (rwebin  rcorein) / nc DTHETA = Pi / (2 * nt) 'Allocate dimension for arrays of radial and polar coordinates ReDim XC((nc + nw + 1) * (nt + 1)), YC((nc + nw + 1) * (nt + 1)), THETAC(nt + 1) 'Allocate dimensions for stiffness arrays ReDim KG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)), PR((nc + nw) * (nt)), PRI((nc + nw) * (nt)) ReDim MG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)), PRTotal((nc + nw) * (nt)), PRStress((nc + nw) * (nt)) 'Gauss points gp(1) = 0.57735 gp(2) = 0.57735 End Sub '********************************************************************************** Sub CartesianMesh() 'generates cartesian mesh 'polar mesh with refining 'geometric mesh refinement in tangential direction stepfactor = 0.92 powerfactor = 0.4 traystep = (stepfactor  1) / (stepfactor ^ nt  1) For q1 = 1 To nt + 1 „THETAC(q1) = THETAC(q1  1) + ((Pi / 2) * (1 / (nt + 1))) THETAC(q1) = (Pi / 2) * (traystep * (stepfactor ^ (q1  1)  1) / (stepfactor  1)) ^ powerfactor Next 'geometric mesh refinement in radial direction ReDim RWEBC(nw + 1) stepfactor = 0.9 powerfactor = 0.8 hstep = (stepfactor  1) / (stepfactor ^ nw  1) For q1 = 1 To nw + 1 RWEBC(q1) = (rwebout  rwebin) * (hstep * (stepfactor ^ (q1  1)  1) / (stepfactor  1)) ^ powerfactor Next 'X,Y coordinates for core For q1 = 1 To nc radpos = rcorein + (q1  1) * hcore For q2 = 1 To (nt + 1) 'nodenumber NN NN = ((nt + 1) * (q1  1)) + q2 XC(NN) = radpos * Cos(THETAC(q2)) YC(NN) = radpos * Sin(THETAC(q2)) Next Next 'X,Y coordinates for web For q1 = nc + 1 To nc + nw + 1 radpos = rwebin + RWEBC(q1  nc) For q2 = 1 To (nt + 1) 'nodenumber NN NN = (nt + 1) * (q1  1) + q2 XC(NN) = radpos * Cos(THETAC(q2)) YC(NN) = radpos * Sin(THETAC(q2)) Next Next Worksheets("FE Nodes").Activate Worksheets("FE Nodes").Cells.ClearContents For q1 = 1 To NN Range(Cells(q1 + 5, 5), Cells(q1 + 5, 5)) = XC(q1) Range(Cells(q1 + 5, 6), Cells(q1 + 5, 6)) = YC(q1) Next End Sub '********************************************************************************** Sub RollWinder() NL = Round((rwebout  rwebin) / hWINDER, 0) 'Read force from input sheet TWI = Range(Cells(23, 2), Cells(23, 2)) TWF = Range(Cells(24, 2), Cells(24, 2)) 'calculate core stiffness for isotropic core CS = Ercore * (rwebin ^ 2  rcorein ^ 2) / (rwebin ^ 2 + rcorein ^ 2  vrtcore * (rwebin ^ 2  rcorein ^ 2)) 'allocate arrays ReDim AW(NL), BW(NL), CW(NL), PW(NL), PIW(NL), tw(NL), RW(NL), ERW(NL) 'mesh generation For q1 = 1 To NL RW(q1) = rwebin + (q1  1) * hWINDER Next 'web line tension profile generation For q1 = 1 To NL tw(q1) = TWI + (q1  1) * (TWF  TWI) / (NL  1) Next 'BEGIN WINDING 'for lap 1 PW(1) = tw(1) * hWINDER / RW(1) ERW(1) = K2 * (PW(1) + K1) 'for lap 2 PW(2) = tw(2) * hWINDER / RW(2) ERW(2) = K2 * (PW(2) + K1) FR = (Etweb / CS  1 + vrtweb) * hWINDER / RW(1) PW(1) = PW(1) + PW(2) / (FR + 1!) ERW(1) = K2 * (PW(1) + K1) 'FINITE DIFFERENCE For q1 = 3 To NL 'lap counter ReDim KW(q1, 3), PIW(q1) For q2 = 2 To q1  1 'form finite difference matrix for layer q2 ERW(q2) = K2 * (PW(q2) + K1) beta = (1  Etweb / ERW(q2)) KW(q2, 3) = (RW(q2) ^ 2 / hWINDER ^ 2 + 3 * RW(q2) / (2 * hWINDER)) KW(q2, 1) = (RW(q2) ^ 2 / hWINDER ^ 2  3 * RW(q2) / (2 * hWINDER)) KW(q2, 2) = (2 * RW(q2) ^ 2 / hWINDER ^ 2 + beta) Next 'boundary conditions for adding lap q1 FR = (Etweb / CS  1 + vrtweb) * hWINDER / RW(1) KW(1, 2) = FR + 1 KW(1, 3) = 1 KW(q1, 2) = 1 'solve for incremental presssures from adding lap q1 'reduction For q2 = 2 To q1 KW(q2, 2) = KW(q2, 2)  KW(q2  1, 3) * KW(q2, 1) / KW(q2  1, 2) Next 'back substiution PIW(q1) = (tw(q1) * hWINDER / RW(q1)) / KW(q1, 2) For q2 = q1  1 To 1 Step 1 PIW(q2) = PIW(q2 + 1) * KW(q2, 3) / KW(q2, 2) Next For q2 = 1 To q1 PW(q2) = PW(q2) + PIW(q2) Next Next '***************COMPRESSION PROBLEM PRESSURE INPUT ReDim QWEBC(nw) For q1 = 1 To nw QWEBC(q1) = Round(NL * (hstep * (stepfactor ^ (q1  1)  1) / (stepfactor  1)) ^ powerfactor, 0) Next For q2 = 1 To nt PR(nc * nt + q2) = PW(Round(QWEBC(2) / 2)) Next For q1 = 2 To nw For q2 = 1 To nt PR(nc * nt + (q1  1) * nt + q2) = PW(QWEBC(q1)) Next Next End Sub '********************************************************************************** Sub StiffnessSystem() 'clear system arrays ReDim KG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)) '*****Dummy loop for filling 0 values (only for checking) For q1 = 1 To 2 * (nc + nw + 1) * (nt + 1) For q2 = 1 To 2 * (nc + nw + 1) * (nt + 1) KG(q1, q2) = 0 Next Next 'form core region stiffness matrix Call PolarCoreMat For q1 = 1 To nc For q2 = 1 To nt x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 Call CartesianMatCon Call StiffnessElement Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 KG(p(m1), p(m2)) = KG(p(m1), p(m2)) + ke(m1, m2) Next Next Next Next 'form web region stiffness matrix For q1 = nc + 1 To nc + nw For q2 = 1 To nt NE = (q1  1) * nt + q2 PRTotal(NE) = PR(NE) Erweb = K2 * (K1  PRTotal(NE)) 'Erweb = Etweb * 0.1 Call PolarWebMat x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 Call CartesianMatCon Call StiffnessElement Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 KG(p(m1), p(m2)) = KG(p(m1), p(m2)) + ke(m1, m2) Next Next Next Next 'Calculate symmetric part of global K matrix 'K matrix is assembled as a square matrix in order to make further calculations easy For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1)) KG(q2, q1) = KG(q1, q2) Next Next '**************'****Add stiffness of cantilevered spring system in rhometer Krod = 15.99 keff = 2 * ((nt + 1) * (nc + nw + 1)) KG(keff, keff) = KG(keff, keff) + Krod '********************************************************************************************************* 'Worksheets("KG").Activate 'Worksheets("KG").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'Range(Cells(5 + q1, 5 + q2), Cells(5 + q1, 5 + q2)) = KG(q1, q2) 'Next 'Next End Sub '********************************************************************************** Sub SystemStiffnessUpdater() 'clear system arrays ReDim KG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)) '*****Dummy loop for filling 0 values (only for checking) 'For q1 = 1 To 2 * (nc + nw + 1) * (nt + 1) 'For q2 = 1 To 2 * (nc + nw + 1) * (nt + 1) 'KG(q1, q2) = 0 'Next 'Next 'form core region stiffness matrix Call PolarCoreMat For q1 = 1 To nc For q2 = 1 To nt x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 Call CartesianMatCon Call StiffnessElement Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 KG(p(m1), p(m2)) = KG(p(m1), p(m2)) + ke(m1, m2) Next Next Next Next 'form web region stiffness matrix For q1 = nc + 1 To nc + nw For q2 = 1 To nt NE = (q1  1) * nt + q2 PRTotal(NE) = PR(NE) + PRStress(NE) Erweb = K2 * (K1  PRTotal(NE)) 'Erweb = Etweb * 0.1 Call PolarWebMat x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 Call CartesianMatCon Call StiffnessElement Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 KG(p(m1), p(m2)) = KG(p(m1), p(m2)) + ke(m1, m2) Next Next Next Next 'Calculate symmetric part of global K matrix 'K matrix is assembled as a square matrix in order to make further calculations easy For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1)) KG(q2, q1) = KG(q1, q2) Next Next '**************Add stiffness of cantilevered spring system in rhometer Krod = 15.99 keff = 2 * ((nt + 1) * (nc + nw + 1)) KG(keff, keff) = KG(keff, keff) + Krod '************************************ End Sub '********************************************************************************** Sub PolarCoreMat() Erase D matcon = Etcore  Ercore * vrtcore ^ 2 D(1, 1) = Ercore * Etcore / matcon D(1, 2) = Ercore * Etcore * vrtcore / matcon D(2, 1) = D(1, 2) D(2, 2) = Etcore ^ 2 / matcon D(3, 3) = Ercore * Etcore / ((Ercore + Etcore) * (1 + vrtcore)) End Sub Sub CartesianMatCon() Erase RD RD(1, 1) = D(1, 1) * Cos(mpa) ^ 4 + 2 * (D(1, 2) + 2 * D(3, 3)) * Sin(mpa) ^ 2 * Cos(mpa) ^ 2 + D(2, 2) * Sin(mpa) ^ 4 RD(1, 2) = (D(1, 1) + D(2, 2)  4 * D(3, 3)) * Sin(mpa) ^ 2 * Cos(mpa) ^ 2 + D(1, 2) * (Cos(mpa) ^ 4 + Sin(mpa) ^ 4) RD(1, 3) = (D(1, 1)  D(1, 2)  2 * D(3, 3)) * Sin(mpa) * Cos(mpa) ^ 3 + (D(1, 2)  D(2, 2) + 2 * D(3, 3)) * Cos(mpa) * Sin(mpa) ^ 3 RD(2, 1) = RD(1, 2) RD(2, 2) = D(2, 2) * Cos(mpa) ^ 4 + 2 * (D(1, 2) + 2 * D(3, 3)) * Sin(mpa) ^ 2 * Cos(mpa) ^ 2 + D(1, 1) * Sin(mpa) ^ 4 RD(2, 3) = (D(1, 1)  D(1, 2)  2 * D(3, 3)) * Sin(mpa) ^ 3 * Cos(mpa) + (D(1, 2)  D(2, 2) + 2 * D(3, 3)) * Sin(mpa) * Cos(mpa) ^ 3 RD(3, 1) = RD(1, 3) RD(3, 2) = RD(2, 3) RD(3, 3) = (D(1, 1) + D(2, 2)  2 * D(1, 2)  2 * D(3, 3)) * Sin(mpa) ^ 2 * Cos(mpa) ^ 2 + D(3, 3) * (Cos(mpa) ^ 4 + Sin(mpa) ^ 4) End Sub '********************************************************************************** Sub StiffnessElement() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 Erase ke For z1 = 1 To 2 For z2 = 1 To 2 t = gp(z1) s = gp(z2) 'Derivatives of shape functions sf1t = (1 + s) / 4 sf1s = (1 + t) / 4 sf2t = (1  s) / 4 sf2s = (1  t) / 4 sf3t = (1  s) / 4 sf3s = (1  t) / 4 sf4t = (1 + s) / 4 sf4s = (1 + t) / 4 'Calculate Jacobian J11 = x1 * sf1t + x2 * sf2t + x3 * sf3t + x4 * sf4t J12 = y1 * sf1t + y2 * sf2t + y3 * sf3t + y4 * sf4t J21 = x1 * sf1s + x2 * sf2s + x3 * sf3s + x4 * sf4s J22 = y1 * sf1s + y2 * sf2s + y3 * sf3s + y4 * sf4s detJ = J11 * J22  J12 * J21 'linear calculate strain gradient matrix BL(1, 1) = (J12 * sf1s + J22 * sf1t) / detJ BL(3, 1) = (J11 * sf1s  J21 * sf1t) / detJ BL(2, 2) = BL(3, 1) BL(3, 2) = BL(1, 1) BL(1, 3) = (J12 * sf2s + J22 * sf2t) / detJ BL(3, 3) = (J11 * sf2s  J21 * sf2t) / detJ BL(2, 4) = BL(3, 3) BL(3, 4) = BL(1, 3) BL(1, 5) = (J12 * sf3s + J22 * sf3t) / detJ BL(3, 5) = (J11 * sf3s  J21 * sf3t) / detJ BL(2, 6) = BL(3, 5) BL(3, 6) = BL(1, 5) BL(1, 7) = (J12 * sf4s + J22 * sf4t) / detJ BL(3, 7) = (J11 * sf4s  J21 * sf4t) / detJ BL(2, 8) = BL(3, 7) BL(3, 8) = BL(1, 7) For m1 = 1 To 8 For m2 = m1 To 8 For m3 = 1 To 3 For m4 = 1 To 3 ke(m1, m2) = ke(m1, m2) + BL(m3, m1) * RD(m3, m4) * BL(m4, m2) * te * detJ Next Next Next Next Next Next End Sub '********************************************************************************** Sub StiffnessConstraint() '*********Apply constraints to Global Stiffness matrix by Elimination approach 'Create dummy stiffness matrix for elimination of constrained DOF ReDim Kdummy(2 * (nw + nc + 1) * (nt + 1), 2 * (nw + nc + 1) * (nt + 1)) For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1) Kdummy(q1, q2) = KG(q1, q2) Next Next 'Identify dof for applying geometric constraints For q3 = 1 To (nc + nw + 1) 'Horizontal axes constraints Yconstraints = 2 * ((nt + 1) * (q3  1) + 1) 'Vertical axes constraints Xconstraints = 2 * ((nt + 1) * (q3))  1 '*********Eliminate constrained dof (delete corresponding rows and columns) XCons = Xconstraints  2 * (q3  1) YCons = Yconstraints  2 * (q3  1) For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) 'Column Eliminations If q2 < YCons And q2 < XCons Then Kdummy(q1, q2) = Kdummy(q1, q2) ElseIf q2 >= YCons And q2 < XCons Then If q2 + 1 = XCons Then Kdummy(q1, q2) = Kdummy(q1, q2 + 2) Else Kdummy(q1, q2) = Kdummy(q1, q2 + 1) End If ElseIf q2 >= XCons Then Kdummy(q1, q2) = Kdummy(q1, q2 + 2) End If Next Next For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) 'Row Eliminations If q1 < YCons And q1 < XCons Then Kdummy(q1, q2) = Kdummy(q1, q2) ElseIf q1 >= YCons And q1 < XCons Then If q1 + 1 = XCons Then Kdummy(q1, q2) = Kdummy(q1 + 2, q2) Else Kdummy(q1, q2) = Kdummy(q1 + 1, q2) End If ElseIf q1 >= XCons Then Kdummy(q1, q2) = Kdummy(q1 + 2, q2) End If Next Next Next '********Calculate constrained global stiffness matrix ReDim Kconstrained(((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))), (2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To ((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q2 = 1 To ((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) Kconstrained(q1, q2) = Kdummy(q1, q2) Next Next 'Worksheets("Kconstrained").Activate 'Worksheets("Kconstrained").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'Range(Cells(5 + q1, 5 + q2), Cells(5 + q1, 5 + q2)) = Kconstrained(q1, q2) 'Next 'Next End Sub '********************************************************************************** Sub MassConstraint() '*********Apply constraints to Global Lumped Mass matrix by Elimination approach 'Create dummy Lumped Mass matrix for elimination of constrained DOF ReDim Mdummy(2 * (nw + nc + 1) * (nt + 1), 2 * (nw + nc + 1) * (nt + 1)) For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1) Mdummy(q1, q2) = MG(q1, q2) Next Next 'Identify dof for applying geometric constraints For q3 = 1 To (nc + nw + 1) 'Horizontal axes constraints Yconstraints = 2 * ((nt + 1) * (q3  1) + 1) 'Vertical axes constraints Xconstraints = 2 * ((nt + 1) * (q3))  1 '*********Eliminate constrained dof (delete corresponding rows and columns) XCons = Xconstraints  2 * (q3  1) YCons = Yconstraints  2 * (q3  1) For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) 'Column Eliminations If q2 < YCons And q2 < XCons Then Mdummy(q1, q2) = Mdummy(q1, q2) ElseIf q2 >= YCons And q2 < XCons Then If q2 + 1 = XCons Then Mdummy(q1, q2) = Mdummy(q1, q2 + 2) Else Mdummy(q1, q2) = Mdummy(q1, q2 + 1) End If ElseIf q2 >= XCons Then Mdummy(q1, q2) = Mdummy(q1, q2 + 2) End If Next Next For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) 'Row Eliminations If q1 < YCons And q1 < XCons Then Mdummy(q1, q2) = Mdummy(q1, q2) ElseIf q1 >= YCons And q1 < XCons Then If q1 + 1 = XCons Then Mdummy(q1, q2) = Mdummy(q1 + 2, q2) Else Mdummy(q1, q2) = Mdummy(q1 + 1, q2) End If ElseIf q1 >= XCons Then Mdummy(q1, q2) = Mdummy(q1 + 2, q2) End If Next Next Next '********Calculate constrained global stiffness matrix ReDim Mconstrained(((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))), (2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To ((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q2 = 1 To ((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) Mconstrained(q1, q2) = Mdummy(q1, q2) Next Next 'Worksheets("Mconstrained").Activate 'Worksheets("Mconstrained").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'Range(Cells(5 + q1, 5 + q2), Cells(5 + q1, 5 + q2)) = Mconstrained(q1, q2) 'Next 'Next End Sub '********************************************************************************** Sub PolarPointer() 'uses dummies q1dummy,q2dummy 'call with q1,q2 loop p(1) = 2 * ((q1  1) * (nt + 1) + q2)  1 p(2) = 2 * ((q1  1) * (nt + 1) + q2) p(3) = 2 * ((q1  1) * (nt + 1) + q2 + 1)  1 p(4) = 2 * ((q1  1) * (nt + 1) + q2 + 1) p(5) = 2 * ((q1) * (nt + 1) + q2)  1 p(6) = 2 * ((q1) * (nt + 1) + q2) p(7) = 2 * ((q1) * (nt + 1) + q2 + 1)  1 p(8) = 2 * ((q1) * (nt + 1) + q2 + 1) End Sub '********************************************************************************** Sub PolarWebMat() Erase D matcon = Etweb  Erweb * vrtweb ^ 2 D(1, 1) = Erweb * Etweb / matcon D(1, 2) = Erweb * Etweb * vrtweb / matcon D(2, 1) = D(1, 2) D(2, 2) = Etweb ^ 2 / matcon D(3, 3) = 2 * Erweb ' * Etweb / ((Erweb + Etweb) * (1 + vrtweb)) End Sub '********************************************************************************** Sub MassSystemLumped() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 'clear system arrays ReDim MG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)) '*****Dummy loop for filling 0 values (only for checking) For q1 = 1 To 2 * (nc + nw + 1) * (nt + 1) For q2 = 1 To 2 * (nc + nw + 1) * (nt + 1) MG(q1, q2) = 0 Next Next '**************************************** 'form core region lumped mass matrix For q1 = 1 To nc For q2 = 1 To nt x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) Call MassCoreElementLumped Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 MG(p(m1), p(m2)) = MG(p(m1), p(m2)) + Mele(m1, m2) Next Next Next Next 'form web region lumped mass matrix For q1 = nc + 1 To nc + nw For q2 = 1 To nt x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) Call MassWebElementLumped Call PolarPointer 'Assemble Global Lumped Mass matrix (Diagnol matrix) For m1 = 1 To 8 For m2 = m1 To 8 MG(p(m1), p(m2)) = MG(p(m1), p(m2)) + Mele(m1, m2) Next Next Next Next 'Worksheets("MG").Activate 'Worksheets("MG").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'Range(Cells(5 + q1, 5 + q1), Cells(5 + q1, 5 + q1)) = MG(q1, q1) 'Next End Sub '********************************************************************************** Sub MassCoreElementLumped() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 Erase Mele For z1 = 1 To 2 For z2 = 1 To 2 t = gp(z1) s = gp(z2) 'Derivatives of shape functions sf1t = (1 + s) / 4 sf1s = (1 + t) / 4 sf2t = (1  s) / 4 sf2s = (1  t) / 4 sf3t = (1  s) / 4 sf3s = (1  t) / 4 sf4t = (1 + s) / 4 sf4s = (1 + t) / 4 'Calculate Jacobian J11 = x1 * sf1t + x2 * sf2t + x3 * sf3t + x4 * sf4t J12 = y1 * sf1t + y2 * sf2t + y3 * sf3t + y4 * sf4t J21 = x1 * sf1s + x2 * sf2s + x3 * sf3s + x4 * sf4s J22 = y1 * sf1s + y2 * sf2s + y3 * sf3s + y4 * sf4s detJ = J11 * J22  J12 * J21 'Calculate elemntal lumped mass matrix for core region For m1 = 1 To 8 Mele(m1, m1) = (te * detJ * rhocore) / (4 * 386) Next Next Next End Sub '********************************************************************************** Sub MassWebElementLumped() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 Erase Mele For z1 = 1 To 2 For z2 = 1 To 2 t = gp(z1) s = gp(z2) 'Derivatives of shape functions sf1t = (1 + s) / 4 sf1s = (1 + t) / 4 sf2t = (1  s) / 4 sf2s = (1  t) / 4 sf3t = (1  s) / 4 sf3s = (1  t) / 4 sf4t = (1 + s) / 4 sf4s = (1 + t) / 4 'Calculate Jacobian J11 = x1 * sf1t + x2 * sf2t + x3 * sf3t + x4 * sf4t J12 = y1 * sf1t + y2 * sf2t + y3 * sf3t + y4 * sf4t J21 = x1 * sf1s + x2 * sf2s + x3 * sf3s + x4 * sf4s J22 = y1 * sf1s + y2 * sf2s + y3 * sf3s + y4 * sf4s detJ = J11 * J22  J12 * J21 'Calculate elemntal lumped mass matrix for core region For m1 = 1 To 8 Mele(m1, m1) = (te * detJ * rhoweb) / (4 * 386) Next Next Next End Sub '********************************************************************************** Sub MassEffective() '****************initial impact on tip node meff = 2 * ((nt + 1) * (nc + nw + 1)) MG(meff, meff) = MG(meff, meff) + (0.00071456) 'Worksheets("MG").Activate 'Worksheets("MG").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'Range(Cells(5 + q1, 5 + q2), Cells(5 + q1, 5 + q2)) = MG(q1, q2) 'Next 'Next End Sub '********************************************************************************** Sub MassLumpedInverse() ReDim MGINV((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)), (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MGINV(q1, q1) = 1 / Mconstrained(q1, q1) Next 'Worksheets("Minv").Activate 'Worksheets("Minv").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'Range(Cells(5 + q1, 5 + q1), Cells(5 + q1, 5 + q1)) = MGINV(q1, q1) 'Next End Sub '********************************************************************************** Sub MinvF() 'Calculate [MGinv] * [F] ReDim MinF((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MinF(q1) = MinF(q1) + (MGINV(q1, q1) * F0(q1)) Next End Sub '********************************************************************************** Sub MinvMU00() 'Calculate [MinM] * [U00] ReDim MinM((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)), (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MinM(q1, q1) = MinM(q1, q1) + (MGINV(q1, q1) * MG(q1, q1)) Next ReDim MinMU00((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MinMU00(q1) = MinMU00(q1) + (MinM(q1, q1) * U00(q1)) Next End Sub '********************************************************************************** Sub TIBracketU0() 'Calculate 2*[M] for bracketed term in deformation eq ReDim MTI2((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)), (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MTI2(q1, q1) = MTI2(q1, q1) + (Mconstrained(q1, q1) * 2) Next 'Multiply (2*[M]*d0) ReDim M2d0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) M2d0(q1) = M2d0(q1) + (MTI2(q1, q1) * U0(q1)) Next 'Calculate delT^2*[K] for bracketed term in deformation eq Call SystemStiffnessUpdater Call StiffnessConstraint ReDim KT2((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)), (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) KT2(q1, q2) = KT2(q1, q2) + (Kconstrained(q1, q2) * delT ^ 2) Next Next 'Multiply (delT^2*[K]*d0) ReDim KT2d0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) KT2d0(q1) = KT2d0(q1) + (KT2(q1, q2) * U0(q2)) Next Next 'Calculate bracketed term {2[M](delT^2)[K]}d0 in deformation eq ReDim BrTI((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) BrTI(q1) = BrTI(q1) + (M2d0(q1)  KT2d0(q1)) Next 'Multiply bracketed term by M Inverse {[Minv]*(2[M](delT^2)[K]})*d0)} ReDim BrMinv((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) BrMinv(q1) = BrMinv(q1) + (MGINV(q1, q1) * BrTI(q1)) Next End Sub '********************************************************************************** Sub TimestepAbaqus() ' Calculate smallest element dimension in the mesh (Differenece between Y coordinate of tip node and its adjacent node) t1 = ((nt + 1) * (nc + nw)) + (nt + 1) t2 = ((nt + 1) * (nc + nw  1)) + (nt + 1) Lmin = Abs(YC(t1)  YC(t2)) 'Calculate dilation wave speed in terms of Lame's constant and material density lamda = (Etweb * vrtweb) / ((1 + vrtweb) * (1  (2 * vrtweb))) mu = Etweb / (2 * (1 + vrtweb)) Cd = Sqr((lamda + (2 * mu)) / (rhoweb / 386)) 'Calculate time step delT = ((Lmin / Cd)) * (0.707) * MultiFactor Worksheets("FE Nodes").Activate Worksheets("FE Nodes").Cells.ClearContents 'For q1 = 1 To NN 'Range(Cells(q1 + 5, 5), Cells(q1 + 5, 5)) = XC(q1) 'Range(Cells(q1 + 5, 6), Cells(q1 + 5, 6)) = YC(q1) 'Next End Sub '********************************************************************************** Sub TransientDI() cons = delT nmax = 25000 ReDim U00((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), U0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim U1((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), Un((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim V00((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), V0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim V1((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), Vn((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim A00((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), A0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim A1((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), An((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim F0(2 * (nc + nw + 1) * (nt + 1)) ReDim Ag(nmax) '*****Erase worksheet data for previous runs Worksheets("Plot").Activate Worksheets("Plot").Cells.ClearContents 'Worksheets("Deformation").Activate 'Worksheets("Deformation").Cells.ClearContents 'Worksheets("Deceleration").Activate 'Worksheets("Deceleration").Cells.ClearContents 'Input initial vectors For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) U0(q1) = 0 'Displacement vector at t=0 V0(q1) = 0 'Velocity vector at t=0 A0(q1) = 0 'Acceleration vector at t=0 F0(q1) = 0 'External force vector at t=0 Next 'Input initial velocity at t=0 invel = (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) V0(invel) = V0(invel)  (32.31) 'Calculate displacment for time step before t=0 For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) U00(q1) = U00(q1)  (delT * V0(q1)) Next Call TIBracketU0 'Call MinvF 'Call MinvMU00 '******Calculate for nmax time steps For n1 = 1 To nmax ReDim Un((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim An((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) '***Un(q1) = Un(q1) + ((MinF(q1) * delT ^ 2) + (((MinMU0(q1) * 2)  (HU0(q1) * delT ^ 2))  (MinMU00(q1)))) Un(q1) = Un(q1) + ((BrMinv(q1))  (U00(q1))) Next Call Udumstress 'Calculate deceleration from deformations For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'An(q1) = An(q1) + ((HUn(q1))) An(q1) = An(q1) + ((Un(q1)  (2 * U0(q1)) + (U00(q1))) / delT ^ 2) Next tipnode = (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) Ag(n1) = An(tipnode) / 386 '****Print results on worksheets for all nodes 'Worksheets("Plot").Activate 'Range(Cells(5 + n1, 3), Cells(5 + n1, 3)) = delT 'Range(Cells(5 + n1, 4), Cells(5 + n1, 4)) = cons 'Range(Cells(5 + n1, 5), Cells(5 + n1, 5)) = Un(tipnode 'Range(Cells(4 + n1, 7), Cells(4 + n1, 7)) = Ag(n1) 'Range(Cells(4 + n1, 9), Cells(4 + n1, 9)) = Rhos '****Reassign values to vectors for carrying out calculations for next time step For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) U00(q1) = U0(q1) U0(q1) = Un(q1) Next Call CartesianMatUpdater Call TIBracketU0 cons = cons + delT Next 'Calculate hardness from maximum deceleration Rhos = 0 For q1 = 1 To nmax If Ag(q1) > Rhos Then Rhos = Ag(q1) End If Next Rhos = Rhos / 2.6 Worksheets("INPUT").Activate Range(Cells(30, 2), Cells(30, 2)) = Rhos End Sub '********************************************************************************** Sub Udumstress() ReDim Udummy(2 * (nw + nc + 1) * (nt + 1)) '2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1) ReDim Ucopy(2 * (nw + nc + 1) * (nt + 1)) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) Ucopy(q1) = Un(q1) Next 'Identify dof for applying geometric constraints For q3 = 1 To (nc + nw + 1) 'Horizontal axes constraints Yconstraints = 2 * ((nt + 1) * (q3  1) + 1) 'Vertical axes constraints Xconstraints = 2 * ((nt + 1) * (q3))  1 '*********Eliminate constrained dof (delete corresponding rows and columns) XCons = Xconstraints '+ 2 * (q3  1) YCons = Yconstraints '+ 2 * (q3  1) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) + (2 * q3) 'Row addition: Add zero deformation to constrained DOF If q1 < YCons And q1 < XCons Then Udummy(q1) = Ucopy(q1) ElseIf q1 = YCons Then Udummy(q1) = 0 ElseIf q1 > YCons And q1 < XCons Then Udummy(q1) = Ucopy(q1  1) ElseIf q1 = XCons Then Udummy(q1) = 0 ElseIf q1 > XCons Then Udummy(q1) = Ucopy(q1  2) End If Next '*****Update Ucopy vector For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) + (2 * q3) Ucopy(q2) = Udummy(q2) Next Next 'Worksheets("Dummy").Activate 'Worksheets("Dummy").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'Range(Cells(5 + q1, 4), Cells(5 + q1, 4)) = Udummy(q1) 'Next 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'Range(Cells(5 + q1, 2), Cells(5 + q1, 2)) = Un(q1) 'Next End Sub '********************************************************************************** Sub CartesianMatUpdater() ReDim PRStress((nc + nw) * (nt)) For q1 = nc + 1 To nc + nw For q2 = 1 To nt NE = (q1  1) * nt + q2 x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 'mean polar angle Erweb = K2 * (K1  PRTotal(NE)) Call PolarWebMat Call CartesianMatCon Call PolarPointer Call CartesianStress '************add stress increment to initial stress PRStress(NE) = PRStress(NE) + (cstress(1) * Cos(mpa) ^ 2 + cstress(2) * Sin(mpa) ^ 2 + 2 * cstress(3) * Cos(mpa) * Sin(mpa)) / 4 Next Next End Sub '********************************************************************************** Sub CartesianStress() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 Erase cstress 'erase cartesian stress components For z1 = 1 To 2 For z2 = 1 To 2 t = gp(z1) s = gp(z2) 'derivatives of shape functions sf1t = (1 + s) / 4 sf1s = (1 + t) / 4 sf2t = (1  s) / 4 sf2s = (1  t) / 4 sf3t = (1  s) / 4 sf3s = (1  t) / 4 sf4t = (1 + s) / 4 sf4s = (1 + t) / 4 'calculate jakobien J11 = x1 * sf1t + x2 * sf2t + x3 * sf3t + x4 * sf4t J12 = y1 * sf1t + y2 * sf2t + y3 * sf3t + y4 * sf4t J21 = x1 * sf1s + x2 * sf2s + x3 * sf3s + x4 * sf4s J22 = y1 * sf1s + y2 * sf2s + y3 * sf3s + y4 * sf4s detJ = J11 * J22  J12 * J21 'linear calculate strain gradient matrix BL(1, 1) = (J12 * sf1s + J22 * sf1t) / detJ BL(3, 1) = (J11 * sf1s  J21 * sf1t) / detJ BL(2, 2) = BL(3, 1) BL(3, 2) = BL(1, 1) BL(1, 3) = (J12 * sf2s + J22 * sf2t) / detJ BL(3, 3) = (J11 * sf2s  J21 * sf2t) / detJ BL(2, 4) = BL(3, 3) BL(3, 4) = BL(1, 3) BL(1, 5) = (J12 * sf3s + J22 * sf3t) / detJ BL(3, 5) = (J11 * sf3s  J21 * sf3t) / detJ BL(2, 6) = BL(3, 5) BL(3, 6) = BL(1, 5) BL(1, 7) = (J12 * sf4s + J22 * sf4t) / detJ BL(3, 7) = (J11 * sf4s  J21 * sf4t) / detJ BL(2, 8) = BL(3, 7) BL(3, 8) = BL(1, 7) For m1 = 1 To 3 For m2 = 1 To 3 For m3 = 1 To 8 cstress(m1) = cstress(m1) + RD(m1, m2) * BL(m2, m3) * Udummy(p(m3)) Next Next Next Next Next End Sub '********************************************************************************** ADVISER‟S APPROVAL: Dr. J. K. Good VITA Dipesh Dilip Mistry Candidate for the Degree of Master of Science Thesis: CONNECTING WINDING CODES TO QUALITY CONTROL DEVICES. Major Field: Mechanical and Aerospace Engineering Biographical: Personal Data: Born in Vasai, India 28th September 1985, son of Mr. Dilip D. Mistry and Mrs. Darshana D. Mistry. Education: Received Diploma in Mechanical Engineering from Bhausaheb Vartak Polytechnic, Vasai in June 2004. Received Bachelor of Science degree from University of Pune, Pune, in May 2007. Completed the requirements for the Master of Science majoring in Mechanical and Aerospace Engineering at Oklahoma State University, Stillwater, Oklahoma in July 2010. Experience: Worked as a teaching assistant for Advanced Experimental Methods in Design (Jan 2008 to May 2010) and Finite Element Methods (Aug 2009 to Dec 2009) at OSU. Worked as a Graduate Research Assistant at Web handling Research Center (WHRC), OSU (June 2008 to July 2010). ADVISER‟S APPROVAL: Dr. J. K. Good Name: Dipesh Dilip Mistry Date of Degree: July, 2010 Institution: Oklahoma State University Location: Stillwater, Oklahoma Title of Study: CONNECTING WINDING CODES TO QUALITY CONTROL DEVICES. Pages in Study: 62 Candidate for the Degree of Master of Science Major Field: Mechanical and Aerospace Engineering Scope and Method of Study: The goal of this study was to develop an extension to the winding models which can convert the outputs of pressures and stresses to units that can be measured using existing quality control devices, the Rho meter in this case. The Rho meter makes dynamic measurement of the hardness of a wound roll. The peak deceleration of the Rho meter striker is measured and converted into arbitrary “Rho” units, which is a measure of hardness of the wound roll. Initially, experiments were conducted to study the variation of hardness in a wound roll as a function of winding tension and the roll outer radius. Based on the experimental results, a dynamic impact model was developed to simulate the Rho meter. The output obtained from the model was compared to the experimental results. Findings and Conclusions: The model made fair predictions of the roll hardness. The material properties vary in time due to the pressures generated by the impact. The explicit method was used in the dynamic analysis since it facilitates the updating of material properties after each solution in time step. Due to this, large computation times are associated with this model. Different algorithms were used to calculate the time step required for the explicit solution and the algorithm employed in this model shows good convergence. The accuracy of the model output increases with an increase in mesh density. ADVISER‟S APPROVAL: Dr. J. K. Good
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Connecting Winding Codes to Quality Control Devices 
Date  20100701 
Author  Mistry, Dipesh Dilip 
Keywords  dynamic impact model, explicit finite element analysis, quality control devices for webs, rhometer, webs, winding codes 
Document Type  
Full Text Type  Open Access 
Abstract  The goal of this study was to develop an extension to the winding models which can convert the outputs of pressures and stresses to units that can be measured using existing quality control devices, the Rho meter in this case. The Rho meter makes dynamic measurement of the hardness of a wound roll. The peak deceleration of the Rho meter striker is measured and converted into arbitrary "Rho" units, which is a measure of hardness of the wound roll. Initially, experiments were conducted to study the variation of hardness in a wound roll as a function of winding tension and the roll outer radius. Based on the experimental results, a dynamic impact model was developed to simulate the Rho meter. The output obtained from the model was compared to the experimental results. The model made fair predictions of the roll hardness. The material properties vary in time due to the pressures generated by the impact. The explicit method was used in the dynamic analysis since it facilitates the updating of material properties after each solution in time step. Due to this, large computation times are associated with this model. Different algorithms were used to calculate the time step required for the explicit solution and the algorithm employed in this model shows good convergence. The accuracy of the model output increases with an increase in mesh density. 
Note  Thesis 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  CONNECTING WINDING CODES TO QUALITY CONTROL DEVICES By DIPESH DILIP MISTRY Bachelor of Science in Mechanical Engineering University of Pune Maharashtra, India 2007 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE July, 2010 iii CONNECTING WINDING CODES TO QUALITY CONTROL DEVICES Thesis Approved: Dr. J. K. Good Thesis Adviser Dr. R. P. Singh Dr. J. D. Jacob Dr. Mark E. Payton Dean of the Graduate College iv ACKNOWLEDGMENTS First of all, I would like to thank the Almighty for giving me the strength during the course of this research. I wish to express my sincere gratitude to Dr. J. K. Good, my adviser, for providing me an opportunity to work for him. He has been a great inspiration and motivated me all the way during my stay at OSU. I thank him for his wisdom and insight. He has been very considerate and supportive of me all the way. I would also like to thank my thesis committee members, Dr. R. P. Singh and Dr. J. D. Jacob. I wish to express my special appreciation to Mr. Ron Markum for guiding me through the experiments and helping me with anything I needed during this research. I thank my friends and colleagues who have always stood by me. Their opinions, suggestions and encouragement have been very helpful in various parts of my research. Finally, I am grateful to my parents, Mr. Dilip D. Mistry and Mrs. Darshana D. Mistry, and my family members, especially my elder brother, Mr. Bhupesh D. Mistry and my sister in law, Mrs. Upasana B. Mistry, who have provided immense psychological and moral support by constantly encouraging and believing in me throughout my life. Dipesh Dilip Mistry.v TABLE OF CONTENTS Chapter ........................................................................................................................... Page I. INTRODUCTION TO WINDING AND QUALITY CONTROL OF WOUND ROLLS ..... 1 II. LITERATURE REVIEW ................................................................................................ 6 2.1 Winding Models: .................................................................................................... 6 2.2 Wound Roll Contact Models: ................................................................................ 9 2.3 Quality Control Devices: .......................................................................................10 2.3.1 Schmidt Hammer: ......................................................................................... 11 2.3.2 PARO Tester: ................................................................................................ 11 2.3.3 Beloit Rho Meter: ..........................................................................................12 2.3.4 Online Density Measurement: ......................................................................12 2.4 Research Objective: ..............................................................................................13 2.5 Accomplishing the Objective: ...............................................................................14 III. EXPERIMENTAL PROCEDURE ................................................................................15 3.1 Winding Setup: .....................................................................................................15 3.2 Hardness Profile Measurement: .......................................................................... 17 3.3 Core Effect on Roll Hardness: ..............................................................................18 3.4 Hardness Measurements in MD and CMD: .........................................................19 3.5 Web Material Properties: .....................................................................................21vi 3.5.1 Stack Test: .....................................................................................................21 3.5.2 Stretch Test: ................................................................................................. 23 3.6 Properties of Polyester web used in this research (PET 377): ............................ 24 IV. DYNAMIC IMPACT MODEL ..................................................................................... 25 4.1 Model Development: ........................................................................................... 25 4.2 Finite Element Formulation: .............................................................................. 26 4.3 Modeling Rho Meter Parameters: ....................................................................... 32 4.4 Transient Analysis: ..............................................................................................40 4.5 Time Step Calculation: ........................................................................................ 42 V. RESULTS AND DISCUSSION .................................................................................... 45 5.1 Experimental Results: .......................................................................................... 45 5.2 Results from the Model: ...................................................................................... 53 5.3 Comparison of Results from Model with Experimental Results: ........................ 55 VI. CONCLUSION AND FUTURE SCOPE ......................................................................60 Future Scope: .............................................................................................................61 BIBLIOGRAPHY ............................................................................................................... 63 APPENDICES ................................................................................................................... 65vii LIST OF FIGURES Figure ............................................................................................................................. Page Figure 1.1: Quality Control Instruments for Wound Rolls ................................................. 4 Figure 3.1: High Speed Winding Line – Web Handling Research Center, OSU ...............16 Figure 3.2: Hardness Measurement of Wound Roll using the Rho Meter ........................ 17 Figure 3.3: Dog Bone Edge Formation due to Thickness Variation in the Web .............. 20 Figure 3.4: Experimental Setup for Stack Test ................................................................. 22 Figure 3.5: PET 377 Web Roll ........................................................................................... 24 Figure 4.1: Finite Element Formulation of Dynamic Impact Model ................................ 28 Figure 4.2: Internal Components of Rho Meter ............................................................... 33 Figure 4.3: Internal Components of Rho Meter (Top View) ............................................ 33 Figure 4.4: Drill Rods and Hinges Supporting the Moving Plastic Parts ......................... 34 Figure 4.5: Setup for Calculating Stiffness of Plastic Hinges ........................................... 35 Figure 4.6: Initial (A) and Final (B) Position of the Spring Mass System ....................... 36 Figure 4.7: Acceleration of the Rho Meter Striker Output on the Oscilloscope ............... 37 Figure 4.8: Initial Position of the Striker at Rest ............................................................. 39 Figure 4. 9: Finite Central Difference Estimates ..............................................................40 Figure 4. 10 Flowchart for Dynamic Impact Model .......................................................... 44 Figure 5.1: Roll Radius Vs Hardness at Constant Winding Tension ................................ 45viii Figure 5.2: Pressure Profile in Rolls Wound at Constant Tension ................................... 46 Figure 5.3: Comparison of Hardness of the Rolls Wound at Constant Tensions with Peak and Average Pressure in the Wound Roll ......................................................................... 47 Figure 5.4: Winding Tension Profile in the Roll due to Alternating Tensions ................. 49 Figure 5.5: Pressure Profile in Rolls Wound at Alternating Winding Tensions ............... 49 Figure 5.6: Comparison of Hardness Profile in the Roll Wound at Constant (5 PLI) and Alternating Tensions (5 PLI  3 PLI) ................................................................................ 50 Figure 5.7: Comparison of Hardness Profile in the Roll Wound at Constant (6 PLI) and Alternating Tensions (2 PLI  6 PLI) .................................................................................51 Figure 5.8: Hardness Comparison of Rolls in MD and CMD (PET 377) .......................... 52 Figure 5.9: Hardness Comparison of Rolls in MD and CMD (Mitsubishi Polyester) ...... 53 Figure 5. 10: Deceleration at Tip Node obtained by using different Time Steps (Δt) in Transient Analysis ............................................................................................................ 54 Figure 5.11: Comparison of Model Output for Linear and Non Linear Material Characteristics with Experimental data ........................................................................... 55 Figure 5.12: Dynamic Model Output for Different Mesh Size, Material Properties are Not Updated during Dynamic Analysis. .................................................................................. 56 Figure 5.13: Variation in Er calculated from Winding Pressures in an Element due to the Size of the Mesh. ................................................................................................................57 Figure 5.14: Dynamic Model Output for Different Mesh Size, Material Properties are Updated during Dynamic Analysis. .................................................................................. 59 ix LIST OF SYMBOLS AND ABBREVIATIONS WHRC Web Handling Research Center HSWL High Speed Web Line PET Polyester Film MD Machine Direction CMD Cross Machine Direction Er Radial Young‟s modulus Eθ Tangential Young‟s modulus εr Radial Strain εθ Tangential Strain σr Radial Stress σθ Tangential Stress σz Axial Stress νrθ Poisson‟s ratio of stress in θdirection to strain in rdirection r Radius of wound roll rc Outer radius of core x P Normal pressure K1 Material constant K2 Material constant C0, C1,..., Cn Constants in Hakiel‟s polynomial for pressure Tw Web tension t Web thickness ρweb Weight density of web material ρcore Weight density of core material h Web thickness made dimensionless by division by core radius δl Change in length te Thickness of quadrilateral element Elemental material matrix in Cartesian coordinates Elemental material matrix in cylindrical coordinates Transformation matrix θj Orientation for jth node [ke] Elemental stiffness matrix [me] Elemental mass matrix [KG], K Global stiffness matrix xi [MG], M Global mass matrix {F}, F (t) Global load vector [H] System matrix B StrainDisplacement matrix J Jacobian η, ξ Natural coordinates βk, βm Penalty number for stiffness and mass matrix respectively Keff Effective stiffness of drill rods Meff Effective mass of Rho meter striker and other components ωn First angular frequency of Rho meter striker d Displacement Velocity Acceleration Δt Time step for transient analysis Lmin Smallest element dimension in the mesh Cd Dilatational wave speed λ Lame‟s constant pli Pounds per linear inch 1 CHAPTER I INTRODUCTION TO WINDING AND QUALITY CONTROL OF WOUND ROLLS Webs can be defined as a material whose length to width ratio is very high and the thickness is very small as compared to both the length and width of the web. Webs are usually manufactured as thin, continuous and flexible layers of material. Webs compose a wide range of materials such as paper, plastic films, textiles, nonwoven materials, metal foils and composites. The most convenient form to store and transport webs is to wind them into rolls. These rolls are then unwound and rewound into subsequent production processes, which collectively adds value to these webs. Winding is the process by which the web is wound on a core to form a roll. The winding operation is done with machines called as winders. There are several types of winders used in the web handling industry. These are classified based on number of drums used. The most common types are single drum, multiple drum and beltreel winders. Based on the method by which the winding torque is applied, winders are categorized as center winders, surface winders, center winder with idler nips, center/surface winders, multiple nip winders, etc. 2 The process of winding is directly affected by the type of winder, the web material and the winder operating parameters (e.g. winding tension, speed, etc). The winding tension results in internal stresses that are responsible for giving the wound roll some integrity as a structure. The quality of the wound roll depends on the stresses which exist in it. Defects such as internal stress bursts, buckling, tearing, etc may occur in the wound roll if these stresses are not properly controlled. Thus, it is a stiff challenge for the web manufacturer and the winder operator to wind good quality rolls of their high quality products. The stresses in the wound roll are relatable to the roll defects. These stresses can be predicted by using winding models. Many winding models have been developed over the years for analyzing the stresses and pressures developed in the roll during winding. The majority of these winding models solve a second order differential equation in pressure written with respect to radius. This equation has non constant coefficients and must be solved numerically to produce results. Typical results consist of pressure and tangential stress profiles with respect to radius within the roll. Thus, winding models can be used to predict the optimum winding conditions for a given web material. Although, the winding model is a powerful quality control tool in the laboratory, it cannot determine the quality of the wound roll on the production floor. A winder operator must be able to check the wound rolls for defects on the shop floor. The most basic tool used for wound roll testing is “the backtenders friend,” a simple short length of stick used to manually strike the roll surface to judge the hardness across the roll width 3 by the sound and feel of the stick impacting the roll surface. With no means to quantify, the method limits the value of this test as a subjective determination of roll structure. (1) The most popular quality control devices on the shop floor are shown in Figure 1.1: Quality Control Instruments for Wound Rolls ,. These are handheld units such as the Beloit1 Rho Meter, the Schmidt2 Hammer and the PAROtester3. These devices measure the hardness of wound roll, which help in detecting the defects during winding. The Beloit Rho Meter, like the backtender‟s stick, is an impact tester. It measures the peak deceleration of a hammer striking the roll surface. The peak deceleration of the striker is converted to an arbitrary unit called “Rho”. The builtin scale permits quantifying the test results. The Schmidt Hammer (2) is another impact tester measuring the rebound of a plunger striking the roll surface. The magnitude of rebound is recorded on a builtin mechanical scale designated in “R” units. The PAROtest (3) is initiated by launching a spring loaded body against the test surface. The impact and rebound velocities are compared resulting in an instantaneous numerical hardness value designated in “L” units. This unit is related to the coefficient of restitution. 1 Beloit Manhattan Division. 1910 Lane Blvd. Kalamazoo, Michigan 49003. USA 2 Testing Machines, Inc. 2 Fleetwood Court, Ronkonkoma, NY 11779 USA 3 Testing Machines, Inc. 2 Fleetwood Court, Ronkonkoma, NY 11779 USA4 Figure 1.1: Quality Control Instruments for Wound Rolls (3), (2) The winding models output the stresses and pressures in the wound roll in terms of engineering units, such as pressure or stress. The quality control devices used on the shop floor output the roll hardness in arbitrary units of hardness. The arbitrary units of hardness cannot be directly coupled to the output in engineering units produced by the winding models. This study focuses on establishing a connection between the outputs of a winding model and a quality control device, the Rho Meter in this case. The Rho meter was selected due to availability in the labs of the Web Handling Research Center at OSU. Initially, experiments were conducted to study the nature of hardness variation in wound rolls as a function of the current outer radius of a winding roll. The relation between the peak pressure in the wound roll and the hardness of the wound roll was investigated. The dynamic hardness measurement range in terms of roll radius was also investigated. These experimental results were used to verify a model which combines a winding model and a dynamic impact model and can output the hardness of the wound roll in terms of 5 Rho units. This model can be an extension to all existing winding models. The experimental procedures and development of the model are discussed in subsequent chapters. The results obtained from the dynamic model are compared to the experimental results. 6 CHAPTER II LITERATURE REVIEW 2.1 Winding Models: Literature on the analysis of wound roll structure first appeared in 1960s. A wound roll structure can be defined as the web layers wound on a core at appropriate winding tension so that the wound roll has some integrity as a structure. Over the years, many winding models were developed for analyzing the stresses and pressures developed in the roll during winding. The early models appeared as 1D models with radial variation only. These models assumed a linear isotropic material and used stress formulas developed from the elasticity solutions of thick cylinders, in order to calculate the internal stresses of a wound roll. Although the 1D models could theoretically capture the accretive pressures which result from winding web layers onto a roll, they did not consider the material complexity due to anisotropy and radial nonlinearity that was required to produce the accurate results. Pfeiffer (4) investigated the radial nonlinearity of wound roll material. He found that the relation between the radial pressure and strain was exponential. 7 {2. 1} He developed a wound roll model in which the radial modulus was state dependent on pressure. {2. 2} Here K1 and K2 are material specific parameters calculated by fitting a curve between pressure versus strain data obtained from stack compression experiments of web material. Altmann (5) provided an analytical solution by assuming constant radial modulus, although different from circumferential modulus. Yagoda (6) was the first one to accurately treat the core boundary condition in wound rolls. Hakiel‟s work (7) is considered to be one of the most important steps in wound roll modeling. He used Altmann‟s formulation and developed a nonlinear anisotropic model. Hakiel incorporated his own model for the radial nonlinearity. He used a polynomial in pressure instead of Pfeiffer‟s exponential form. {2. 3} Hakiel‟s model was the first one to treat the state dependent properties and all boundary conditions rigorously. He derived a second order differential equation to calculate the increments in pressure ( ) in the wound roll. 8 {2. 4} He required two boundary conditions for the solution and so he utilized the outer and inner boundary conditions. The inner boundary condition is the compatibility of deformation of the core and the first layer of wound roll {2. 5} The outer boundary condition comes from the thin pressure vessel formula {2. 6} Hakiel employed a finite difference scheme to solve the set of equations for the addition of each lap and calculate the incremental pressure in each lap. The incremental pressures in each lap are summed with all previous increments to determine the total pressure. The material properties ( ) due to the current total pressure were then calculated and updated in the model prior to solving for the increments in pressure due to the addition of the next layer. After all layers are added the tangential stresses as a function of radius can be determined using the equilibrium equation posed in cylindrical coordinates: {2. 7} 9 2.2 Wound Roll Contact Models: Wound roll contact models have been developed over the years to analyze the contact problems that occur during winding or storage of wound roll. Ganapathi (8) developed a 2D axisymmetric finite element model to study the diametral compaction of wound roll, with state dependent properties, between rigid surfaces. He calculates the initial pressures in the wound roll by using Hakiel‟s (7) nonlinear winding model. The winding pressure is then used to calculate the radial modulus by using Pfeiffer‟s (4) relation. The wound roll cylinder is treated as an anisotropic cylinder for the formulation of the finite element contact model. The cylinder is meshed by using 4 noded quadrilateral elements (9). He adopted the general strategy of node to node contact in order to find the force which will bring down an associated node to the contacting surface. The contact is accomplished via linear interpolation wherein it requires three solutions per contacting node. The first solution is obtained by applying a unit load at the core and calculating the associated pull down force at the node. The next pull down force is obtained by applying half unit load and the final solution is obtained by calculating the load by linear interpolation of the first two results such that it will produce zero pull down force. At this point, the stresses are calculated and the material properties are updated. As a result of the contact pressures combining with the pressures due to winding, the radial modulus will now vary spatially with respect to radius and tangential location within the contact model. The operation is terminated when the total compression force is obtained. He verified his model experimentally (8). Mollamahmutoglu (10) developed a Nip and Wound Roll Contact model to analyze the contact pressures between nip and wound roll. The model is based on 2D plane stress finite element formulation of the contact of two orthotropic cylinders, one of 10 them represents the wound roll with nonlinear material properties and other represents the rigid nip roller. Similar to Ganapathi‟s model, the initial pressures due to winding are calculated by using Hakiel‟s (7) nonlinear winding model. The winding pressure is then used to calculate the radial modulus by using Pfeiffer‟s (4) relation. Mollamahmutoglu uses a very simple and efficient algorithm based on quasilinearization technique to obtain a rapid solution. The algorithm proposed by Mollamahmutoglu is much faster because it solves only once per contact node whereas the algorithm adopted by Ganapathi‟s model solves three times per node. This model was verified experimentally using Ganapathi‟s data and good comparison was found with Ganapathi‟s model when the problem of contact between wound roll and rigid wall was solved using the same technique. 2.3 Quality Control Devices: Smith (11) contends that roll hardness is probably the most important factor in determining the difference between a good and bad roll. Rolls that are wound too soft may have defects like starring and telescoping due to interlayer slippage during unwinding. Rolls that are wound too tight contain high inwound tension due to which tension bursts can occur inside the wound roll. The main aim of wound roll quality control measurement techniques is to determine if visually indiscernible defects are present into the roll. If the defects are detected then the rolls are rewound and future defects can be prevented by changing the parameters which affect the wound roll structure. Most of the quality control devices are handheld units and measure the roll hardness on the principle of impact. A few quality control 11 devices will be discussed briefly in this chapter. A thorough study of the Rho meter was conducted as we need to understand the instrument for modeling it into the dynamic impact finite element model that we intend to develop as a part of this research. This will be discussed in subsequent chapters. 2.3.1 Schmidt Hammer: The Schmidt test hammer is a handheld device based on the principle of rebound measurement. The Schmidt Hammer is composed of a spring loaded impact plunger which is pressed against the wound roll causing the compression of the spring. After reaching a certain point, the spring is released and an internal hammer mass is launched against the impact plunger. The amount of hammer mass rebound depends on the hardness of the wound roll. The magnitude of rebound is recorded on a builtin mechanical scale designated in “R” units. Wound roll hardness profile is generated by taking hardness measurements across the roll width. The hardness profiles are then used to detect the rolls with inconsistent hardness. In addition to the builtin scale, a paper tape print out is also available which gives an instant picture and records the cross roll hardness profile. (2), (1) 2.3.2 PARO Tester: PARO tester is an impact tester which is initiated by launching a spring loaded body against the wound roll surface. The impact and rebound velocities are measured and converted to output an instantaneous numerical hardness value. The value is the quotient of the impact and rebound velocities multiplied by 1000, results in output designated in “L” units. This unit is related to the coefficient of restitution. The accuracy and ease of implementation makes it one of the popular instruments on the shop floor. The digital display and inherent data memory makes it easy to interpret 12 and operate. The measurement information can be downloaded on a computer for further analysis. (3), (1) 2.3.3 Beloit Rho Meter: The Rho meter is an impact tester which measures the peak deceleration of a striker when it impacts on the roll surface. The striker is composed of steel and is a constant mass. A cantilever spring system with straight line motion accelerates the striker mass to a known constant velocity. A peak force is generated by the contact of striker with the surface of the wound roll. The peak force is detected by a high accuracy accelerometer, which is rigidly mounted on the striker to measure the peak deceleration on impact. The peak deceleration is converted into an arbitrary unit called “Rho”, which is a measure of roll hardness. The Rho meter has a digital display which indicates the roll hardness instantly and accurately. It is easy to operate and interpret. The Rho meter produces highly repeatable readings for a given roll. (4), (1) 2.3.4 Online Density Measurement: David Roisum‟s PhD work (12) was on online density measurement. The paper industry has been making online density measurements on winders for many years. It is easy to measure by using two encoders. One encoder is attached to an external nip and it measures the amount of web entering the wound roll and the peripheral velocity of the wound roll. A second encoder is attached to the core and measures the angular velocity of the core. During a measurement period one can determine 1. What length of web has entered the winder and, 2. How has the radius of the wound roll changed. 13 Given a web of known thickness one can determine its compressed volume and hence its density. Roisum tried to connect online density measurement to wound roll stresses through winding models. What he found at that time (1990) was that with encoder technology of that time he was marginally able to measure density with the accuracy he needed to predict wound roll deformations on materials as compressible as paper. Conceivably it would have worked better yet on materials with even greater compressibility like tissue. The stress condition inside the wound roll is a key factor in determining the quality of the roll. The winding process and associated parameters are generally responsible for the wound roll‟s stress condition. Stresses induced during winding are known to be related to defects. Therefore, determining the winding parameters that produce optimum roll stress profiles is of importance. Winding models can predict the stresses in terms of engineering units, which can be verified experimentally in a laboratory. On the production floor, the winder operator uses handheld quality control devices to detect the defects in the wound roll by measuring the roll hardness. 2.4 Research Objective: Review of literature shows that only one attempt has been previously made to provide outputs from wound roll models that can also be measured in the laboratory or on the production floor. That attempt was only marginally successful. The goal of this research is to develop extensions to wound roll models which converts outputs of pressures and stresses to units that can be measured using existing quality control devices, the Rho Meter in this case. 14 2.5 Accomplishing the Objective: The variation of roll hardness with respect to the roll geometry and winding conditions is investigated experimentally. The relation between the hardness and peak pressures in the wound roll is also studied. Experiments are also conducted to investigate the dynamic measurement range of the Rho Meter. The internal mechanism of the Rho Meter is studied and modeled into a finite element dynamic impact model with nonlinear material properties. The output from the finite element dynamic impact model is compared to the experimental results. The experimental procedure and development of dynamic impact model is discussed in detail in subsequent chapters. 15 CHAPTER III EXPERIMENTAL PROCEDURE In this research the hardness of wound rolls is measured using a Rho meter. Hardness variation of the rolls with respect to the outer radius is investigated. Rolls are wound at different tension levels. The winder is capable of holding the tension at zero velocity, thus facilitating the hardness measurements at intermediate outer roll radius positions from the core to the final wound roll radius. 3.1 Winding Setup: The High Speed Winding Line (HSWL) consists of number of rolls through which the web passes from unwind roller to the winding roller. It can reach high speeds of winding, by which the machine gets its name. The tension is applied using an automatic feedback tension control system. The winding machine is equipped with position guides which track and control the lateral position of the web during winding. The surface speed of the web is maintained low enough to avoid air entrainment problem. The webs were wound at five different levels of tension viz. 2, 3, 4, 5 and 6 PLI tension. 16 Figure 3.1: High Speed Winding Line – Web Handling Research Center, OSU 17 Figure 3.1 shows the HSWL. The rolls are wound on a steel core which has known properties and high stiffness. The procedure to measure the hardness of the roll is discussed in the following section. 3.2 Hardness Profile Measurement: Experiments were conducted on the web supplied by DuPont. The material was polyester 377 (PET 377) film in a 6inch width and 200 gage (0.002 inch) thickness. Since the HSWL is capable of holding tension at zero velocity, about half an inch (pile height) of web was wound on the steel core at a constant winding tension and then the machine was set to run at zero velocity. The Rho meter was used to measure the roll hardness at the current pile height of the web wound onto a core. Figure 3.2: Hardness Measurement of Wound Roll using the Rho Meter 18 Figure 3.2 shows the hardness measurement being taken by a Rho meter on a wound roll mounted on the HSWL. The hardness measurement was taken at five locations along the width of the web. The average hardness of the roll for the current pile height was calculated from these readings. The radius of the steel core was 1.7 inches and the final outer radius of the roll was 5.2 inches. The procedure was repeated for increments of one half inch of pile height. Thus, hardness readings were taken at seven radial locations in the wound roll. The hardness readings were plotted as a function of roll radius. The average hardness of the roll wound at different winding tension is compared to the peak pressure in the roll, which is calculated using a winding model. 3.3 Core Effect on Roll Hardness: The Rho meter makes dynamic measurement of hardness over a radial range. Experiments were conducted in order to investigate the range over which the hardness measurements are made by the Rho meter. The methodology adopted to study this range is by winding the web on a steel core (rigid) and then winding it on a comparatively soft core. The comparison of hardness at corresponding radial locations will give us an insight about the range of the Rho meter. We define “Core effect” as the variation in hardness of a roll due to winding it either on a steel (rigid) core or a softer core, at the same winding tension. The web was wound on a steel core at two different tension levels. These tensions were switched after about half inch increment in the pile height. The initial pile (0.5 inch), wound at low tension, acts as a soft core for the preceding pile (0.5 inch) of web that is wound at a higher tension. Experiments were conducted for two cases. For the 19 first case the higher and lower tension levels were 5 PLI and 3 PLI, respectively. For the second case, the tension levels were changed to 6 PLI for higher tension and 2 PLI for lower tension. In any case, as compared to the rigid steel core, each increment in pile height acts as a core for the preceding increment in pile to be wound onto the roll. Initially, the web was wound at high tension starting at the core and then the tension was reduced after one half inch pile height was wound on. The tension was alternated after every one half inch of pile height. After winding one half inch (pile height) of web, the HSWL is held at zero velocity while holding tension. The Rho meter was used to measure the roll hardness across the roll width and the average roll hardness was calculated at the current pile height of the web. Next, winding was started at low tension at the core and then the tension was increased after one half inch pile height. The tension was then alternated after every half inch of pile height. After winding one half inch (pile height) of web, the HSWL is held at zero velocity while holding tension. The Rho meter was used to measure the roll hardness across the roll width and the average roll hardness was calculated at the current pile height of the web. The final outer radius of the roll was 5.2 inches and hardness readings were taken at seven radial locations in the wound roll. The hardness readings were plotted as a function of roll radius. The hardness profile is compared to the pressure profile of the roll, generated in the winding model. 3.4 Hardness Measurements in MD and CMD: The striking mass of the Rho meter has a finite length. It is critical to determine the effect of contact length of the striker on the roll hardness measurements, when it 20 impacts the wound roll. The materials used for these experiments were the PET 377 previously described and a Mitsubishi polyester film in a 24inch width and 300 gage (0.003 inch) thickness. Mitsubishi intentionally produced this web with a highly radical thickness profile which results into very high variation in hardness across the width of the wound roll. Figure 3.3 shows the dog bone edge formation due to thickness variation in a Mitsubishi polyester wound roll. Figure 3.3: Dog Bone Edge Formation due to Thickness Variation in the Web The PET 377 web was wound on the HSWL at three different levels of tension viz. 1.5, 2 and 4 PLI. After winding about two inches (pile height) of web, the HSWL is held at zero velocity while holding tension. Now hardness measurements were taken by holding the Rho meter in the machine direction (MD) and the cross machine direction (CMD) at three locations along the width of the web. The Mitsubishi polyester web was 21 wound on HSWL at a tension of 1.5 PLI. After winding about ten inches (pile height) of web, the HSWL is held at zero velocity while holding tension. Again hardness measurements were taken by holding the Rho meter in MD and CMD at five locations along the width of the web. The hardness data obtained for both cases (MD and CMD) was plotted as a function of roll width. 3.5 Web Material Properties: 3.5.1 Stack Test: The stack test is conducted to measure the radial modulus (Er) of the web material. The radial modulus is a very critical variable in the winding process. Parameters like air entrainment, surface roughness and coatings combine to form the radial modulus of the wound roll, which is state dependent on the pressure in the wound roll. This nonconstant material property offers a challenge in understanding and modeling the wound roll structure. Figure 3.4 shows the experimental setup for stack test to measure the radial modulus (Er) of the web. 22 Figure 3.4: Experimental Setup for Stack Test The web is compressed slowly at constant velocity to allow any air to escape laterally form the stack as the compression progresses. Deformation and compression load are simultaneously recorded using a data acquisition system, during these tests. This data is used to infer the pressure as a function of strain in the stack. The experimental data is then used to produce a pressure versus strain chart. A least square approximation is used to fit a curve to the experimental data and the calculate K1 and K2 in Pfeiffer‟s exponential expression {2.1} for the pressure versus strain curve. K1 has the units of pressure and K2 is dimensionless. 23 3.5.2 Stretch Test: The stretch test is conducted to measure the tangential modulus (Eθ) of the web material. Tangential modulus of the web is an important factor in determining the pressures developed inside the wound roll. Both Er and Eθ affect the radial stresses and thus the pressure in the roll calculated using the second order differential equation {2.4}. The tangential modulus for any common engineering specimen could be measured using a material testing machine. For web materials better results are obtained by stretching a longitudinal length of the web (50 feet) and recording the load required to achieve a measured deformation. This method is preferred because the web must somehow be gripped to install a tension. Grip effects typically induce two dimensional stresses in the vicinity of the grips and thus the stress state in the web is not uniaxial. The long specimen minimizes the errors due to grip effects. In the stretch tests that are done 50 feet of web that is 6 inches wide are stretched. This results in an aspect ratio of 100, which serves to minimize the variability due to grip effects and specimen misalignment. One end of the web is fixed by taping it to the floor and a known load is applied at the other end. The load is constantly increased and the change in length ( δl ) is measured for application of each load. The loading is done for about twenty increments. The strain can then be calculated from the change in length for each load. The stress can be calculated from the area of the web and the force applied. The stress  strain curve is generated and the slope of this curve gives us the tangential modulus (Eθ). 24 3.6 Properties of Polyester web used in this research (PET 377): The 200 gage films have the following measured properties: Web width = 6 inches Web thickness = 0.002 inches Material constants for Pfeiffer‟s (4) relation calculated experimentally: K1 = 3.19 psi K2 = 38.35 Eθ = Ez = 711,000 psi Figure 3.5: PET 377 Web Roll 25 CHAPTER IV DYNAMIC IMPACT MODEL 4.1 Model Development: In this chapter a model is developed for a 2Dplane formulation of the dynamic impact of Rho meter striker and the wound roll. The model consists of two main components; a Winding model and a Dynamic Impact model. The winding model is a development similar to that of Hakiel (7) and the impact model is based upon a dynamic finite element formulation. Both models are developed using Excel and Visual Basic for Applications and the Excel/VBA code is shown in appendix A. The winding model employed in this analysis is based on Hakiel‟s non linear model for wound roll stresses (7). The input needed for this model are material properties which include the tangential modulus ( ) of the web, Young‟s modulus for the core, Pfeiffer‟s constants K1 and K2 (4) to establish the radial modulus of the web, Poisson‟s ratio for the web material and the core and the weight density of the web material and the core. The other input parameters include the geometry of the core and the roll, which include the inner and outer radii of the core, the outer radius of the roll, the thickness of the web and the width of the web material. The winding tension variation with radius is also required. 26 Winding formulations such as that proposed by Hakiel (7) solve for the increments in radial pressure within each layer of the roll due to the addition of the most recent layer. This is done by solving a second order differential equation written in terms of the incremental pressures which relies on two boundary conditions. One boundary condition is that the outside of the core and the inside of the roll must maintain continuity while the second condition assumes that the pressure beneath the outermost layer is known from the winding tension input and the current radius of the outer lap. Once the incremental pressures have been solved for the addition of the most recent layer they are summed with the incremental pressures within the laps due to the addition of all the previous layers to produce the total pressure within each lap. The total pressure in each lap is then used to update the radial direction modulus ( ) of the web using Pfeiffer‟s expression: , where P is the pressure due to winding. With updated properties a new layer is added and the second order differential equation is solved for the addition of the next layer until all layers have been added. After the final lap is added the final winding pressure and the radial modulus are known as functions of wound radius and will be used as initial values for the dynamic impact finite element model. 4.2 Finite Element Formulation: Figure 4.1 shows the finite element model for Rho meter striker impact with a wound roll. The finite element model is simplified by taking advantage of the symmetry and thus only a quarter of the roll is modeled. This simplification is justified because the striker impact effects are localized in the impact region only, which will be verified later, and due to the symmetry offered by the geometry of the wound roll. In this model each 27 node will be allowed to deform in the X and Y directions (u and v respectively), except for those nodes which may be constrained on the X and Y axes as shown in Figure 4.1. 2D quadrilateral plane stress elements with anisotropic properties are used in the formulation. Due to the fact that most of the deformation will occur near the impact region, a mesh refinement technique will be used along the radial and tangential direction that will densify the mesh in the impact region. The finite element formulation begins with mesh generation. The mesh generation consists of two parts. In the initial part a uniform mesh is generated in the radial direction for the core structure and in the second part a refined mesh is generated in the radial direction for the wound roll region. A refined mesh is generated in the tangential direction for the core and the wound roll region (10). 28 Figure 4.1: Finite Element Formulation of Dynamic Impact Model The initial wound roll material properties were computed by the winding model which determined Er as a function of r and E was assumed constant. These properties are now needed in the global Cartesian coordinates of the impact model. Thus the finite element formulation requires transformation of element material properties that are measured in the cylindrical coordinate (rθ) system of the winding model to the Cartesian coordinate (XY) system used in the finite element model. This transformation 29 is employed separately for each finite element during the calculation for the material matrix [ ]. {4. 1} Here is the orthogonal matrix which represents the transformation (10). For the element occupying a (i, j) position it can be given as: {4. 2} The stiffness matrix for each element [ke] is computed by using a 2x2 Gaussian quadrature scheme to numerically integrate (9): {4. 3} where J is the Jacobian and has the form: {4. 4} and, {4. 5} 30 {4. 6} {4. 7} {4. 8} These elemental stiffness matrices are assembled into a Global Stiffness Matrix [KG] by using the direct assembly procedure. The Global Stiffness Matrix [KG] is always symmetric and banded matrix. Although [KG] stored in a banded form saves some computer memory, in this case [KG] is stored as a square matrix in order to reduce the complexity of computations involved in the dynamic analysis. A lumped mass formulation was employed to compute the elemental mass matrix [me] for each element. The total element mass in each direction is distributed equally to the nodes of the element, and the masses are associated with translational degrees of freedom only (9). For the 2D quadrilateral element, the lumped element mass matrix [me] is: {4. 9} Further, the elemental mass matrices are assembled into a Global Lumped Mass Matrix [MG] by using the direct assembly procedure. The mesh refinement in the impact 31 region is very high and thus the lumped mass formulation is fairly accurate in this region. By using the lumped mass formulation a diagonal Global Lumped Mass Matrix [MG] is obtained, which offers computational efficiency in the dynamic analysis algorithm that will be discussed later. Geometric constraints must be enforced to the finite element model, which occur due to symmetry. As seen in figure 4.1, we need to constraint the vertical movement of the nodes on the lower surface (Y=0) and the horizontal movement on the left boundary (X=0) of the quarter cylinder. Initially, the penalty method was used to enforce these geometric boundary conditions. A penalty number (βk) is added to the appropriate components of the global stiffness matrix [KG] and an inertial penalty number (βm) is added to the appropriate components of the global mass matrix [MG]. The penalty approach resulted in huge numbers on the major diagonal of the global stiffness matrix, which required the value of critical time step (Δt) for the transient analysis to be reduced accordingly. Thus the overall time period over which the deceleration of the striker mass is calculated, reduced significantly and the solution never reached the peak deceleration. Due to the difficulties encountered by the penalty approach, the elimination approach was employed for enforcing the geometric constraints. In this method, the constraints are enforced by eliminating the rows and columns in the global stiffness matrix [KG] and the global mass matrix [MG], corresponding to the constrained degree of freedom. After enforcing the constraints, we obtain a reduced global stiffness and mass matrix. The advantage of the elimination approach over the penalty approach is that the value of critical time step (Δt) is increased resulting in a longer time step for transient analysis and decreased computational time. The Rho meter dynamic elements must also be included into the finite element model which is discussed in the following section. 32 4.3 Modeling Rho Meter Parameters: The internal mechanical components of a Rho meter are shown in figures 4.2 and 4.3 below. The striker impacts the wound roll at a known velocity and accelerometer measures the peak deceleration of the striker after the impact. The striker has a straight line motion which is accelerated by a cantilevered spring system to a desired contact velocity. The most important parameters that affect the peak deceleration on impact are the effective stiffness of the cantilevered spring system, the effective mass of the system, the initial velocity of striker and the hardness of the wound roll. For simplicity, the entire Rho meter is not modeled. The important Rho meter parameters are included at the tip node location as shown in Figure 4.1. Due to this, there are no external forces acting on the system and thus the global force vector {FG} is always equal to zero through time in the dynamic model. 33 Figure 4.2: Internal Components of Rho Meter Figure 4.3: Internal Components of Rho Meter (Top View) 34 The effective stiffness of the cantilevered spring system comes from the stiffness of the two steel drill rods and the stiffness of the composite hinges inserted in the moving plastic parts. These components are shown in Figure 4.4. The two drill rods measure 4.75" in length and 0.095" in diameter. These rods are made from steel. With known dimensions and material properties, the stiffness of each of the two drill rods was calculated from equation 4.10; {4. 10} Figure 4.4: Drill Rods and Hinges Supporting the Moving Plastic Parts 35 The stiffness of the hinges supporting the moving plastic parts was determined experimentally by measuring the displacement of the spring mass system and the corresponding force by using a handheld force gage. Figure 4.5 shows the setup for measuring the force and the corresponding displacement of the spring mass system. Figure 4.5: Setup for Calculating Stiffness of Plastic Hinges 36 The stiffness of the drill rods come into consideration after the spring mass system is displaced by 0.0625 inches. This is a nonstandard configuration for the Rho meter due to some rubber grommets that were missing in one of the WHRC‟s two units. The force applied by a handheld force gage up to a displacement of 0.0625 inches is directly affected by the stiffness of the hinges supporting the moving plastic mass. Figure 4.6 shows the initial and final positions of the spring mass system before and after the force is applied by the force gage. Figure 4.6: Initial (A) and Final (B) Position of the Spring Mass System The force required for producing a deformation of 0.0625 inches was found to be 0.58 lb. The stiffness of the hinges was calculated by using these readings; {4. 11} Thus, the effective stiffness of the system will be the sum of the stiffnesses of the drill rods and the plastic hinges. 37 {4. 12} The striker of the Rho meter is expected to stay in contact with the wound roll surface during the deceleration period. Thus, the effective stiffness ( ) of the system is simply added to the tip mode location of the global stiffness matrix [KG]. The natural angular frequency ( ) of the system can be measured experimentally by connecting the output of the accelerometer to an oscilloscope. The acceleration sensed by the accelerometer is output on the oscilloscope to measure the natural frequency of the system. The natural frequency of the system is calculated from the two consecutive peak acceleration observed after the impact of the striker. Figure 4.7 shows the acceleration of the Rho meter striker output on the oscilloscope. Figure 4.7: Acceleration of the Rho Meter Striker Output on the Oscilloscope 38 The natural angular frequency ( ) of the system observed on the oscilloscope was 23.80 Hz. Thus; {4. 13} The effective mass ( ) of the system is a summation of the mass of the striker, the support rod, the aluminum spring support, the moving plastic mass, the accelerometer, and the effective mass of the cantilevered drill rods shown in figures 4.2 and 4.3. The effective mass ( ) of the system can be calculated from equation; {4. 14} The effective mass ( ) of the system is simply added to the tip mode location of the global lumped mass matrix [MG]. The initial velocity of the Rho meter striker can be calculated from the first natural angular frequency ( ) of the system and assuming a harmonic displacement function. A harmonic displacement function (d) was assumed of the form: {4. 15} {4. 16} when the striker is released, ; where do is the height at which the striker is released above the wound roll by the Rho meter instrument. When is , the velocity will be maximum and will be zero. 39 {4. 17} Figure 4.8 shows the initial position of the striker at rest. As observed in Figure 4.8, the striker is not in line with the surface of the Rho meter shoe plate at time (t) =0. The distance by which the striker needs to displace in order to strike the surface of the wound roll was measured to be 0.125 inches. The time required for this displacement of the striker was calculated to be 7.016 x 103 sec. Figure 4.8: Initial Position of the Striker at Rest Thus the velocity of the striker on impact with the surface of the wound roll will be; 40 {4. 18} The velocity of the striker on impact ( ) is added to the tip node location of the global velocity vector at the start of the transient analysis. With the Rho meter parameters calculated and input within the appropriate locations in the global stiffness matrix, global mass matrix and global velocity vector, the transient analysis can be performed which is discussed in the following section. 4.4 Transient Analysis: The Direct Integration method is used to perform the transient analysis of the impact of Rho meter striker and the wound roll. In this method the derivatives of the deformation through time are evaluated by using Finite Central Difference method. The central difference estimates of the velocity and acceleration vectors at time ti are calculated using the following relations; Figure 4. 9: Finite Central Difference Estimates 41 {4. 19} {4. 20} The calculation of the time step for the transient analysis is very critical because it directly affects the accuracy of the solution. The calculation of time step will be discussed in the subsequent section. In the dynamic code, the peak deceleration is sought after the impact. The procedure applied in the dynamic code for calculating the decelerations is as follows: 1. Calculate time step 2. Read initial conditions: , which are known. 3. Solve for previous time step; {4. 21} 4. With determined, now solve for ; (for next ) {4. 22} This equation results from inserting 4.20 into the equation of motion shown in 4.22, (13): {4. 23} 5. With given and determined, we now solve for {4. 24} 6. Calculate corresponding deceleration using {4.15} or : 42 {4. 25} 7. Use steps 56 in a recursive fashion until the first peak deceleration is obtained. Once the peak deceleration of the striker is detected by using the above procedure, it is converted into hardness in terms of “Rho‟s”. The conversion is obtained by using the following relation that was used by the manufacturer of the Rho meter: {4. 26} The result, in terms of Rho hardness, is then printed on the Excel sheet. 4.5 Time Step Calculation: The explicit dynamic finite element analysis is very sensitive to the selection of the time step . The explicit procedure integrates through time by using small time increments as discussed in the preceding section. The value of Δt is very critical, because it affects the results directly. If too high a value for Δt is entered, divergent output results. If too low a value for Δt is entered, large rounding errors can result. The time step is generally expressed in terms of the highest natural frequency of the system as given in the following equation; {4. 27} The highest natural frequency ( ) can be calculated by using vector iteration method. In this method, if we use the system matrix , then the solution is converges to the highest natural frequency. The iteration process requires huge 43 computational times for the calculation of time step. Thus we switched to a more efficient method discussed below. A second alternative method allows the time step to be approximated at the elemental level. The time step ( ) for the dynamic problem is calculated based upon the smallest transit time of a dilatational wave across any of the elements in the mesh (14). This results in saving the computational time associated with the vector iteration method. The approximate value for can be calculated from equation {4.28}; {4. 28} Where, {4. 29} „Lmin‟ is the smallest element dimension in the mesh and „Cd‟ is the dilatational wave speed expressed in terms of Lame‟s constant and material density. This estimate for Δt is only approximate and is typically scaled by a factor between 1 and for two dimensional analysis and between 1 and for three dimensional analysis (14). For this case the value obtained by the method discussed above needs to be further scaled by a factor between 0 and 1 in order to compensate for the material non linearity parameters. The value of this factor is dependent on the mesh density of the wound roll. A flow chart of the dynamic impact finite element model is shown in figure 4.10.44 If Dtip=1st Peak Dec YES A Apply geometrical constraints Begin Transient Analysis. 1) Calculate Δt 2) Read Initial Velocity 3) Form initial deformation vector Use Finite – Central Difference Method to estimate the deformations. Calculate acceleration from deformations at previous time steps Calculate deformations at tip node for next Δt Calculate corresponding peak deceleration (Dtip) at tip node Calculate elemental stresses and corresponding radial pressures. Update Global Stiffness matrix. Apply Rho meter parameters & Constraints Calculate Rho Hardness Print Output =Rho Hardness NO STOP START Input: Material Properties, winding parameters, core and roll geometry. Execute Hakiel type Winding model and calculate Pr & Er for the wound roll. Generate Finite element mesh, using mesh refinement in the web region Generate elemental stiffness matrices using Er due to winding. Assemble Global stiffness matrix. Generate elemental mass matrices. Assemble Global Lumped mass matrix. Model Rho meter parameters: Add effective stiffness of the cantilevered spring system to tip node location in the Global Stiffness matrix. Add effective mass of striker to tip node location in the Global Mass matrix. A Figure 4. 10 Flowchart for Dynamic Impact Model 45 CHAPTER V RESULTS AND DISCUSSION 5.1 Experimental Results: The PET 377 web was wound at different winding tensions viz. 2 PLI to 6 PLI. The Rho meter was used to measure the hardness of the roll across the width. The average roll hardness was calculated at the current pile height (radius) of the web. The measurements were taken at seven radial locations. Figure 5.1 shows the hardness profile for the rolls wound at different tensions. Figure 5.1: Roll Radius Vs Hardness at Constant Winding Tension 0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 6 Hardness (Rho) Roll Radius (inch) 2 PLI 3 PLI 4 PLI 5 PLI 6 PLI Winding Tension46 The hardness values obtained at lower winding tensions is low and it increases as the winding tension is increased. The hardness near the core is higher as compared to the hardness away from the core, although, this trend becomes less pronounced at higher winding tensions. It was observed that at low winding tensions (23 PLI), the difference in hardness between consecutive winding tensions is greater as compared to the difference in hardness of rolls wound at higher tensions (46 PLI). Thus it can be inferred that there is a nonlinear relation between the wound roll hardness and the winding tensions. The hardness profile was compared to the pressure profile of the wound roll. Figure 5.2 shows the pressure profile for the rolls wound at different tensions. Figure 5.2: Pressure Profile in Rolls Wound at Constant Tension 0 50 100 150 200 250 300 350 400 450 500 0 1 2 3 4 5 6 Presssure (psi) Roll Radius (Inch) 2 PLI 3 PLI 4 PLI 5 PLI 6 PLI Winding Tension47 The pressures inside the wound roll were generated by using Hakiel‟s winding model. The pressure profile obtained from the winding model is of a plateaued roll where the pressures may be nearly constant throughout the roll. It was observed that, unlike the hardness, the pressure increases exponentially with an increase in winding tension. Figure 5.3: Comparison of Hardness of the Rolls Wound at Constant Tensions with Peak and Average Pressure in the Wound Roll Figure 5.3 shows comparison of the hardness of rolls wound at constant tension and the peak and average pressures in the roll. The pressure in the roll is generated by using Hakiel‟s winding model. The peak pressure occurs at the core of the wound roll, whereas the average pressure is calculated from the pressures at different radial locations within the wound roll. It was observed that the peak pressure increases 0 10 20 30 40 50 60 70 80 0 50 100 150 200 250 300 350 400 450 500 0 1 2 3 4 5 6 7 Hardness (Rho) Pressure (psi) Winding Tension (PLI) Peak Pressure Average Pressure Hardness: Experiments Roll Outer Radius = 5.2 inch48 exponentially with an increase in winding tension. The average pressure in a wound roll is comparatively lower than the peak pressure, but shows a similar trend with an increase in winding tension. The roll hardness shows a different trend. The slope of the hardness curve decreases with an increase in winding tension. Thus it can inferred that after a certain tension level, the hardness measurements obtained by the Rho meter will reach an asymptote. As observed in Figure 5.1, the hardness readings of the roll near the core are higher due to the effect of the core. The core in this case is radially much stiffer than the web layers in compression. This core effect is observed because the Rho meter makes dynamic measurement of the roll hardness over a localized range. Rolls were wound by alternating the winding tension after every one half inch pile height and the hardness measurements are taken across the roll width. Figure 5.4 shows the winding tension profile in the rolls wound at alternating tension. The pressure profile for these rolls was generated using Hakiel‟s winding model and is shown Figure 5.5. Experiments were conducted for two cases. For the first case the higher and lower tension levels were 5 PLI and 3 PLI, respectively. For the second case, the tension levels were changed to 6 PLI for higher tension and 2 PLI for lower tension. 49 Figure 5.4: Winding Tension Profile in the Roll due to Alternating Tensions Figure 5.5: Pressure Profile in Rolls Wound at Alternating Winding Tensions 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 Winding Tension (PLI) Roll Radius (inch) Alternating Tensions 5 PLI 3 PLI Alternating Tensions 2 PLI 6 PLI STEEL CORE Web Roll 0 30 60 90 120 150 180 210 240 270 300 330 360 0 1 2 3 4 5 6 Pressure (psi) Roll Radius (Inch) Alternating Tension 5 PLI 3 PLI Alternating Tension 2 PLI 6 PLI50 The hardness profiles for the rolls wound at alternating tension are shown in Figure 5.6 and Figure 5.7. It was observed that the average hardness of the roll wound on a steel core is higher when compared to the average hardness of the rolls wound on a core, whose stiffness is lower than the steel core. Thus, it can be inferred that the Rho meter hardness measurements are sensitive to changes in winding tension. With the data available from this experiment, it is difficult to find the radial range over which the Rho meter makes the hardness measurements. In order to investigate the “Core Effect” and the radial range of the Rho meter, the hardness measurements need to be made at relatively lower pile heights than half of an inch. These measurements would be effective in capturing the actual trend of the roll hardness profile in the radial direction. The experimental results can be used to study the variation in hardness of the wound roll due to stiffness of the core. Figure 5.6: Comparison of Hardness Profile in the Roll Wound at Constant (5 PLI) and Alternating Tensions (5 PLI  3 PLI) 0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 6 Rho Hardness (Rho) Roll Radius (inch) Alternating Tensions 5 PLI 3 PLI 5 PLI Constant Tension51 Figure 5.7: Comparison of Hardness Profile in the Roll Wound at Constant (6 PLI) and Alternating Tensions (2 PLI  6 PLI) Experiments were conducted to determine the effect of contact length of the Rho meter striker on the roll hardness measurements, when it impacts the wound roll. Two different web materials were used for this experiment, PET 377 and Mitsubishi polyester. PET 377 web was wound at different winding tensions and hardness readings were taken across the width of the web for each winding tension. The hardness measurements were taken by holding the Rho meter in the machine direction (MD) and the cross machine direction (CMD) at these locations. 0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 6 Rho Hardness (Rho) Roll Radius (inch) Alternating Tensions 2 PLI 6 PLI 6 PLI Constant Tension52 Figure 5.8: Hardness Comparison of Rolls in MD and CMD (PET 377) Figure 5.8 and 5.9 show the comparison of hardness of the roll for measurements taken in MD and CMD at the same locations along the width. In Figure 5.9 it can be observed that the hardness of the wound roll is very high at one of the edges. This results from the formation of dog bone edge due to very high thickness variation in the web. 0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 Hardness (Rho) Roll width (inch) 1.5 PLI MD 2 PLI MD 4 PLI MD 1.5 PLI CMD 4 PLI CMD 2 PLI CMD Winding Tension Roll OuterRadius = 4.2 inches53 Figure 5.9: Hardness Comparison of Rolls in MD and CMD (Mitsubishi Polyester) Good comparison was observed for the hardness measurements taken in MD and CMD. Thus it can be inferred that the contact length of the striker on impact has a less pronounced effect on the hardness measurements. 5.2 Results from the Model: The winding model and the finite element dynamic model is developed using Excel and Visual Basic for Applications. Initially, the convergence of the model was verified by comparing the outputs obtained for a constant set of input data by varying the time step (Δt) for transient analysis. Three different values for Δt were used for the transient analysis and the deceleration at the tip node was plotted as a function of time 0 20 40 60 80 100 120 140 160 0 5 10 15 20 25 Hardness (Rho) Roll width (inch) 1.5 PLI MD 1.5 PLI CMD Winding Tension Roll OuterRadius = 10 inches54 for all the three cases. Figure 5. 10 shows the deceleration at the tip node obtained by using three different values for Δt in the transient analysis. The deceleration curve from the model has a similar trend and the peak deceleration occurs at the same time for all the three cases. Thus it can concluded that the model shows good convergence when the given time step (Δt) is less than the critical time step for the problem. Figure 5. 10: Deceleration at Tip Node obtained by using different Time Steps (Δt) in Transient Analysis 55 5.3 Comparison of Results from Model with Experimental Results: The model output was compared to the experimental data. The non linear material properties were employed into the model and the material properties were updated after each solution in time step. Figure 5.11: Comparison of Model Output for Linear and Non Linear Material Characteristics with Experimental data Figure 5.11 shows the comparison of model output with the experimental data for two different rolls wound at constant tensions of 2Pli and 6 Pli respectively. The hardness output in terms of Rho units is obtained from the model at different outer roll radius locations. The output obtained from the model for lower winding tension (2 Pli) is 0 10 20 30 40 50 60 70 0 1 2 3 4 5 6 Hardness (Rho) Roll Radius (inch) Experiment 2 PLI Non Linear (Model) 2 PLI Experiment 6 PLI Non Linear(Model) 6 PLI56 in better comparison with the experimental data as compared to the model output for higher winding tension (6 Pli). This difference may be due to the coarse mesh size used for the above simulations. During the dynamic analysis the material properties are updated after each solution in time step, which results in huge computational times. The coarse mesh size was used in order to reduce the computational time for dynamic analysis. Figure 5.12: Dynamic Model Output for Different Mesh Size, Material Properties are Not Updated during Dynamic Analysis. 57 Figure 5.12 shows the model output for different mesh sizes. In order to reduce the computation times, the material properties are not updated during the dynamic analysis. As the mesh size increases, the deceleration output from the model decreases. Typically in FEA increased mesh size leads to results that are more accurate but in fact this is not true in this case. This is due to the estimation of the radial modulus (Er) of the elements in the wound roll region is calculated from the average pressures obtained from the winding model. Figure 5.13: Variation in Er calculated from Winding Pressures in an Element due to the Size of the Mesh. Figure 5.13 shows the variation in Er, calculated from the winding pressures, due to the size of the element in the mesh. As the mesh size gets coarser, the size of the element in the mesh gets larger and thus the average pressure within the element increases resulting in a higher value of Er for the element. Thus, the stiffness of the elements in the coarse mesh is higher as compared to the stiffness of the elements in the fine mesh. The 58 higher stiffness of the elements in the web region results in higher values for deceleration of the Rho meter striker during the dynamic analysis. Thus, if the material properties are not updated during the dynamic analysis, the output obtained from the model for different mesh sizes will not be accurate and huge errors result due to neglecting the change in material properties due to localized increase in pressures during the impact between the Rho meter striker and the wound roll. Next, the output from the dynamic model was obtained for different mesh sizes by updating the material properties after each solution in time step during the dynamic analysis. Due to this the computation times were high for this model and thus the deceleration output from the model was obtained only for a few time steps. Figure 5.14 shows the model output for different mesh sizes. As the mesh size increases, the deceleration output from the model increases. This will produce model results in Figure 5.11 that are closer to the test results. This is due to the fact that the Er for the elements in the web region is a function of the total pressure (PrTotal) due to winding and the localized pressures due to impact of the Rho meter striker. 59 Figure 5.14: Dynamic Model Output for Different Mesh Size, Material Properties are Updated during Dynamic Analysis. As the element size becomes smaller, the average pressures calculated at the four Gauss points within the element will be more accurate. The total pressure (PrTotal) in the element will result into higher Er for the element and thus the stiffness of the element will increase with a finer mesh size. These results support the argument that a finer mesh sizes will produce more accurate results from the model. 60 CHAPTER VI CONCLUSION AND FUTURE SCOPE In this research the experimental data was utilized to develop a model for dynamic impact of the Rho meter striker with the surface of the wound roll. From the results obtained from the experiments we can conclude that; The hardness variation in the wound roll with respect to the winding tension and the outer roll radius is studied experimentally. The hardness of the wound roll increases with an increase in the winding tension. Pressures and Hardness increase with an increase in tension, but the slope of the curves show a different trend. It can be inferred that the hardness measurements obtained by the Rho meter will increase up to a certain point, with an increase in winding tension, after that the measurements obtained by the Rho meter will reach an asymptote. The Rho meter makes a dynamic measurement of roll hardness over a range of wound roll radius close to the outer surface where the impact occurred. The output from the device is sensitive to variations in winding tension and the stiffness of the core. Further experimentation is required to study what portion of the outside of the winding roll affect the measurements that are made by the Rho meter. 61 The hardness measurements made by the Rho meter by aligning the axis of the Rho meter striker in MD and CMD are similar. Thus the geometry of impact of the Rho meter striker against the wound roll does not affect the hardness measurements by a considerable amount. The geometric constraints enforced by using the Elimination approach works better than the geometric constraints enforced by using the Penalty approach for this case. Although, the elimination approach increases the complexity due to change in dimension of the related matrices, they result into lower computation times. The algorithm used in the model for calculating the time step (Δt) shows good convergence for a given set of inputs. The selection of Δt is influenced greatly by the mesh density. The assumptions made in the development of the dynamic impact model are reasonable as the model results are in fair agreement with the experimental results. The output obtained from the dynamic model is in better comparison with the experimental data only if the material properties are updated locally during the dynamic analysis. Although the time required for computations is very high in this case the properties must be updated to produce accurate results. The output from the model obtained by employing a fine mesh size is more accurate as compared to the output obtained by employing a coarse mesh size. Future Scope: The main idea behind this study was to develop an extension to the winding models which can convert the outputs of pressures and stresses to units that can be measured using existing quality control devices. We have successfully modeled the Rho 62 meter in this research but it has been shown that reasonable mesh densities and thus large computation times are required to produce good converged results. To this end the computational times of the model can be reduced by using a localized impact model instead of using a quarter model of the wound roll, as modeled in this case. The region affected by the impact can be determined by examining the variation in pressures of the elements. The region where the pressure variation due to impact is negligible can be omitted and only the affected region could be modeled. This will result in more accurate results with reduced computational times for the dynamic analysis. To get more confidence in this model, experiments similar to the ones we conducted in this study need to be conducted for different web material and the results should be compared with the output from the model. This work could be extended further to develop similar models for some other popular quality control devices like the Schmidt hammer and the PARO tester. 63 BIBLIOGRAPHY 1. Wound Roll Testers. Paper Industry Web (PIW). [Online] www.paperindustryweb.com/rolltest.htm. 2. Paper Roll Hardness Tester. Schmidt Roll Hardness Tester. [Online] Testing Machines Inc. http://www.testingmachines.com/4020paperrollhardnesstesterrecording.html. 3. Paper Roll Hardness Tester. PARO Tester2. [Online] Testing Machine Inc. http://www.testingmachines.com/4022parotester2.html. 4. Pfeiffer, J. D., Internal Pressures in a Wound Roll of Paper. 5, August 1966, TAPPI Journal, Vol. 49, pp. 342  347. 5. Altmann, H. C., Formulas for Computing the Stresses in CenterWound Rolls. 4, April 1968, TAPPI Journal, Vol. 51, pp. 176  179. 6. Yagoda, H. P., Resolution of Core Problem in Wound Rolls. 4, Transactions of ASME, December 1980, Journal of Applied Mechanics, Vol. 47, pp. 847  854. 7. Hakiel, Z., Nonlinear Model for Wound Roll Stresses. 5, May 1987, TAPPI Journal, vol 70, no. 5, Vol. 70, pp. 113  117. 8. Ganapathi, S., Diametral Compression of Wound Rolls with State Dependent Properties. MS Thesis, Oklahoma State University, December 2007.64 9. Chandrupatla, T. R. and Belegundu, A. D., Introduction to Finite Elements in Engineering. 3. Prentice Hall, Inc, 2003. 10. Mollamahmutoglu, C. A 2D AxisSymmetric Wound Roll Model Including Nip Effects. PhD Dissertation, Oklahoma State University, December 2009. 11. Smith, R. D., The Art of Winding Quality Rolls. August 2001, Paper Film & Foil Converter, pp. 46  53. 12. Roisum, D. R., The Measurement of Roll Stresses During Winding. PhD Dissertation, Oklahoma State University, 1990. 13. Komzsik, L., Computational Techniques of Finite Element Analysis. CRC Press, 2005. pp. 195  202. 14. Kandadai, B., The Development of WoundOnTension in Webs Wound into Rolls. PhD Dissertation, Oklahoma State University, December 2006. APPENDICES APPENDIX A The dynamic impact model developed using Excel and Visual Basic for Applications is shown below: '**********************************************************************************'‟*********DYNAMIC IMPACT MODEL FOR CONNECTING WINDING CODES******** ‟*************************TO QUALITY CONTROL DEVICES************************ '**********************************************************************************' Option Explicit '**************INPUT PARAMETERS Dim nw, nc, nt As Integer Dim rcorein, rwebin, rwebout, te As Double Dim Erweb, Etweb, vrtweb, Ercore, Etcore, vrtcore, rhoweb, rhocore As Double Dim K1, K2 As Double Dim alpha, beta1, MultiFactor As Double Dim start '**************VARIABLES FOR MESH REFINEMENT AND OTHER CALCULATIONS Dim stepfactor, powerfactor, traystep, hstep, QWEBC(), RWEBC() As Double Dim DTHETA, hweb, hcore, Pi, chmax, chmin As Double Dim PR(), PRI(), um(), PRTotal(), PRStress() As Double Dim XC(), YC(), radpos As Double Dim THETAC() As Double Dim p(8) As Integer Dim gp(2) As Double Dim NN As Integer Dim KG(), FG(), Kdummy(), Kconstrained(), Mdummy(), Mconstrained() As Double Dim mpa, x1, x2, x3, x4, y1, y2, y3, y4, RD(3, 3) As Double Dim xi, xo, yi, yo, ri, ro, THETAE As Double Dim J11, J12, J21, J22, BL(3, 8), matcon As Double Dim ke(8, 8), D(3, 3) As DoubleDim t, s, detJ As Double Dim sf1t, sf2t, sf3t, sf4t, sf1s, sf2s, sf3s, sf4s As Double Dim q1, q2, q3, m1, m2, m3, m4, z1, z2, n1, t1, t2, t3 As Integer Dim MG(), CG(), MGINV() As Double Dim Mele(8, 8), NTN(8, 8) As Double Dim PC As Double Dim mpc, meff, keff, Xconstraints, Yconstraints, XCons, YCons As Integer Dim Irod, Drod, Erod, Lrod, Krod As Double '**************WINDER VARIABLES Dim KW(), AW(), BW(), CW(), PW(), PIW(), tw(), RW(), ERW() As Double Dim beta, FR, TWI, TWF, CS, hWINDER As Double Dim NL As Integer '**************DYNAMIC PROBLEM VARIABLES Dim delT, cons, Rhos As Double Dim U00(), U0(), U1(), Un(), Udummy(), Ucopy(), H() As Double Dim V00(), V0(), V1(), Vn() As Double Dim A00(), A0(), A1(), An(), Ag() As Double Dim F0(), MinF(), MinM(), KUC(), G() As Double Dim invel, tipnode, vec1, nmax As Integer Dim NE As Integer Dim cstress(3) As Double Dim HU0(), HUn(), MinMU0(), MinMU00() As Double Dim C, Omega, Iter1(), Iter2(), delTmax, delTnl As Double Dim Lmin, Cd, lamda, mu As Double Dim MTI2(), HmDelT(), BrTI(), BrTDI(), BrMinv(), M2d0(), KT2d0(), KT2() As Double '**************CHARTS 'Dim LoaddefChart, PressureChart, WPressureChart, NCZChart, MPCChart, NodeChart As Object Dim length As Double '********************************************************************************** '********************************************************************************** Sub Main() start = Timer Call Inputdata Application.ScreenUpdating = False Application.DisplayAlerts = False Call CartesianMesh Call RollWinder Call StiffnessSystem Call MassSystemLumped Call MassEffective Call StiffnessConstraint Call MassConstraint Call MassLumpedInverse Call TimestepAbaqus Call TransientDI 'Call CartesianResults Worksheets("INPUT").Activate Range(Cells(31, 2), Cells(31, 2)) = Timer  start End Sub '********************************************************************************** Sub Inputdata() 'mesh parameters and geometrical data hWINDER = Range(Cells(13, 2), Cells(13, 2)) nw = Range(Cells(4, 2), Cells(4, 2)) nt = Range(Cells(5, 2), Cells(5, 2)) nc = Range(Cells(6, 2), Cells(6, 2)) rcorein = Range(Cells(8, 2), Cells(8, 2)) rwebin = Range(Cells(9, 2), Cells(9, 2)) rwebout = Range(Cells(10, 2), Cells(10, 2)) te = Range(Cells(11, 2), Cells(11, 2)) 'material properties Etweb = Range(Cells(16, 2), Cells(16, 2)) vrtweb = Range(Cells(17, 2), Cells(17, 2)) Ercore = Range(Cells(19, 2), Cells(19, 2)) Etcore = Range(Cells(20, 2), Cells(20, 2)) vrtcore = Range(Cells(21, 2), Cells(21, 2)) rhocore = Range(Cells(26, 2), Cells(26, 2)) rhoweb = Range(Cells(27, 2), Cells(27, 2)) MultiFactor = Range(Cells(28, 2), Cells(28, 2)) K1 = Range(Cells(14, 2), Cells(14, 2)) K2 = Range(Cells(15, 2), Cells(15, 2)) 'Pi constant Pi = 3.141592654 'Calculate layer thickness for core and web regions hweb = (rwebout  rwebin) / nw hcore = (rwebin  rcorein) / nc DTHETA = Pi / (2 * nt) 'Allocate dimension for arrays of radial and polar coordinates ReDim XC((nc + nw + 1) * (nt + 1)), YC((nc + nw + 1) * (nt + 1)), THETAC(nt + 1) 'Allocate dimensions for stiffness arrays ReDim KG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)), PR((nc + nw) * (nt)), PRI((nc + nw) * (nt)) ReDim MG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)), PRTotal((nc + nw) * (nt)), PRStress((nc + nw) * (nt)) 'Gauss points gp(1) = 0.57735 gp(2) = 0.57735 End Sub '********************************************************************************** Sub CartesianMesh() 'generates cartesian mesh 'polar mesh with refining 'geometric mesh refinement in tangential direction stepfactor = 0.92 powerfactor = 0.4 traystep = (stepfactor  1) / (stepfactor ^ nt  1) For q1 = 1 To nt + 1 „THETAC(q1) = THETAC(q1  1) + ((Pi / 2) * (1 / (nt + 1))) THETAC(q1) = (Pi / 2) * (traystep * (stepfactor ^ (q1  1)  1) / (stepfactor  1)) ^ powerfactor Next 'geometric mesh refinement in radial direction ReDim RWEBC(nw + 1) stepfactor = 0.9 powerfactor = 0.8 hstep = (stepfactor  1) / (stepfactor ^ nw  1) For q1 = 1 To nw + 1 RWEBC(q1) = (rwebout  rwebin) * (hstep * (stepfactor ^ (q1  1)  1) / (stepfactor  1)) ^ powerfactor Next 'X,Y coordinates for core For q1 = 1 To nc radpos = rcorein + (q1  1) * hcore For q2 = 1 To (nt + 1) 'nodenumber NN NN = ((nt + 1) * (q1  1)) + q2 XC(NN) = radpos * Cos(THETAC(q2)) YC(NN) = radpos * Sin(THETAC(q2)) Next Next 'X,Y coordinates for web For q1 = nc + 1 To nc + nw + 1 radpos = rwebin + RWEBC(q1  nc) For q2 = 1 To (nt + 1) 'nodenumber NN NN = (nt + 1) * (q1  1) + q2 XC(NN) = radpos * Cos(THETAC(q2)) YC(NN) = radpos * Sin(THETAC(q2)) Next Next Worksheets("FE Nodes").Activate Worksheets("FE Nodes").Cells.ClearContents For q1 = 1 To NN Range(Cells(q1 + 5, 5), Cells(q1 + 5, 5)) = XC(q1) Range(Cells(q1 + 5, 6), Cells(q1 + 5, 6)) = YC(q1) Next End Sub '********************************************************************************** Sub RollWinder() NL = Round((rwebout  rwebin) / hWINDER, 0) 'Read force from input sheet TWI = Range(Cells(23, 2), Cells(23, 2)) TWF = Range(Cells(24, 2), Cells(24, 2)) 'calculate core stiffness for isotropic core CS = Ercore * (rwebin ^ 2  rcorein ^ 2) / (rwebin ^ 2 + rcorein ^ 2  vrtcore * (rwebin ^ 2  rcorein ^ 2)) 'allocate arrays ReDim AW(NL), BW(NL), CW(NL), PW(NL), PIW(NL), tw(NL), RW(NL), ERW(NL) 'mesh generation For q1 = 1 To NL RW(q1) = rwebin + (q1  1) * hWINDER Next 'web line tension profile generation For q1 = 1 To NL tw(q1) = TWI + (q1  1) * (TWF  TWI) / (NL  1) Next 'BEGIN WINDING 'for lap 1 PW(1) = tw(1) * hWINDER / RW(1) ERW(1) = K2 * (PW(1) + K1) 'for lap 2 PW(2) = tw(2) * hWINDER / RW(2) ERW(2) = K2 * (PW(2) + K1) FR = (Etweb / CS  1 + vrtweb) * hWINDER / RW(1) PW(1) = PW(1) + PW(2) / (FR + 1!) ERW(1) = K2 * (PW(1) + K1) 'FINITE DIFFERENCE For q1 = 3 To NL 'lap counter ReDim KW(q1, 3), PIW(q1) For q2 = 2 To q1  1 'form finite difference matrix for layer q2 ERW(q2) = K2 * (PW(q2) + K1) beta = (1  Etweb / ERW(q2)) KW(q2, 3) = (RW(q2) ^ 2 / hWINDER ^ 2 + 3 * RW(q2) / (2 * hWINDER)) KW(q2, 1) = (RW(q2) ^ 2 / hWINDER ^ 2  3 * RW(q2) / (2 * hWINDER)) KW(q2, 2) = (2 * RW(q2) ^ 2 / hWINDER ^ 2 + beta) Next 'boundary conditions for adding lap q1 FR = (Etweb / CS  1 + vrtweb) * hWINDER / RW(1) KW(1, 2) = FR + 1 KW(1, 3) = 1 KW(q1, 2) = 1 'solve for incremental presssures from adding lap q1 'reduction For q2 = 2 To q1 KW(q2, 2) = KW(q2, 2)  KW(q2  1, 3) * KW(q2, 1) / KW(q2  1, 2) Next 'back substiution PIW(q1) = (tw(q1) * hWINDER / RW(q1)) / KW(q1, 2) For q2 = q1  1 To 1 Step 1 PIW(q2) = PIW(q2 + 1) * KW(q2, 3) / KW(q2, 2) Next For q2 = 1 To q1 PW(q2) = PW(q2) + PIW(q2) Next Next '***************COMPRESSION PROBLEM PRESSURE INPUT ReDim QWEBC(nw) For q1 = 1 To nw QWEBC(q1) = Round(NL * (hstep * (stepfactor ^ (q1  1)  1) / (stepfactor  1)) ^ powerfactor, 0) Next For q2 = 1 To nt PR(nc * nt + q2) = PW(Round(QWEBC(2) / 2)) Next For q1 = 2 To nw For q2 = 1 To nt PR(nc * nt + (q1  1) * nt + q2) = PW(QWEBC(q1)) Next Next End Sub '********************************************************************************** Sub StiffnessSystem() 'clear system arrays ReDim KG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)) '*****Dummy loop for filling 0 values (only for checking) For q1 = 1 To 2 * (nc + nw + 1) * (nt + 1) For q2 = 1 To 2 * (nc + nw + 1) * (nt + 1) KG(q1, q2) = 0 Next Next 'form core region stiffness matrix Call PolarCoreMat For q1 = 1 To nc For q2 = 1 To nt x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 Call CartesianMatCon Call StiffnessElement Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 KG(p(m1), p(m2)) = KG(p(m1), p(m2)) + ke(m1, m2) Next Next Next Next 'form web region stiffness matrix For q1 = nc + 1 To nc + nw For q2 = 1 To nt NE = (q1  1) * nt + q2 PRTotal(NE) = PR(NE) Erweb = K2 * (K1  PRTotal(NE)) 'Erweb = Etweb * 0.1 Call PolarWebMat x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 Call CartesianMatCon Call StiffnessElement Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 KG(p(m1), p(m2)) = KG(p(m1), p(m2)) + ke(m1, m2) Next Next Next Next 'Calculate symmetric part of global K matrix 'K matrix is assembled as a square matrix in order to make further calculations easy For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1)) KG(q2, q1) = KG(q1, q2) Next Next '**************'****Add stiffness of cantilevered spring system in rhometer Krod = 15.99 keff = 2 * ((nt + 1) * (nc + nw + 1)) KG(keff, keff) = KG(keff, keff) + Krod '********************************************************************************************************* 'Worksheets("KG").Activate 'Worksheets("KG").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'Range(Cells(5 + q1, 5 + q2), Cells(5 + q1, 5 + q2)) = KG(q1, q2) 'Next 'Next End Sub '********************************************************************************** Sub SystemStiffnessUpdater() 'clear system arrays ReDim KG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)) '*****Dummy loop for filling 0 values (only for checking) 'For q1 = 1 To 2 * (nc + nw + 1) * (nt + 1) 'For q2 = 1 To 2 * (nc + nw + 1) * (nt + 1) 'KG(q1, q2) = 0 'Next 'Next 'form core region stiffness matrix Call PolarCoreMat For q1 = 1 To nc For q2 = 1 To nt x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 Call CartesianMatCon Call StiffnessElement Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 KG(p(m1), p(m2)) = KG(p(m1), p(m2)) + ke(m1, m2) Next Next Next Next 'form web region stiffness matrix For q1 = nc + 1 To nc + nw For q2 = 1 To nt NE = (q1  1) * nt + q2 PRTotal(NE) = PR(NE) + PRStress(NE) Erweb = K2 * (K1  PRTotal(NE)) 'Erweb = Etweb * 0.1 Call PolarWebMat x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 Call CartesianMatCon Call StiffnessElement Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 KG(p(m1), p(m2)) = KG(p(m1), p(m2)) + ke(m1, m2) Next Next Next Next 'Calculate symmetric part of global K matrix 'K matrix is assembled as a square matrix in order to make further calculations easy For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1)) KG(q2, q1) = KG(q1, q2) Next Next '**************Add stiffness of cantilevered spring system in rhometer Krod = 15.99 keff = 2 * ((nt + 1) * (nc + nw + 1)) KG(keff, keff) = KG(keff, keff) + Krod '************************************ End Sub '********************************************************************************** Sub PolarCoreMat() Erase D matcon = Etcore  Ercore * vrtcore ^ 2 D(1, 1) = Ercore * Etcore / matcon D(1, 2) = Ercore * Etcore * vrtcore / matcon D(2, 1) = D(1, 2) D(2, 2) = Etcore ^ 2 / matcon D(3, 3) = Ercore * Etcore / ((Ercore + Etcore) * (1 + vrtcore)) End Sub Sub CartesianMatCon() Erase RD RD(1, 1) = D(1, 1) * Cos(mpa) ^ 4 + 2 * (D(1, 2) + 2 * D(3, 3)) * Sin(mpa) ^ 2 * Cos(mpa) ^ 2 + D(2, 2) * Sin(mpa) ^ 4 RD(1, 2) = (D(1, 1) + D(2, 2)  4 * D(3, 3)) * Sin(mpa) ^ 2 * Cos(mpa) ^ 2 + D(1, 2) * (Cos(mpa) ^ 4 + Sin(mpa) ^ 4) RD(1, 3) = (D(1, 1)  D(1, 2)  2 * D(3, 3)) * Sin(mpa) * Cos(mpa) ^ 3 + (D(1, 2)  D(2, 2) + 2 * D(3, 3)) * Cos(mpa) * Sin(mpa) ^ 3 RD(2, 1) = RD(1, 2) RD(2, 2) = D(2, 2) * Cos(mpa) ^ 4 + 2 * (D(1, 2) + 2 * D(3, 3)) * Sin(mpa) ^ 2 * Cos(mpa) ^ 2 + D(1, 1) * Sin(mpa) ^ 4 RD(2, 3) = (D(1, 1)  D(1, 2)  2 * D(3, 3)) * Sin(mpa) ^ 3 * Cos(mpa) + (D(1, 2)  D(2, 2) + 2 * D(3, 3)) * Sin(mpa) * Cos(mpa) ^ 3 RD(3, 1) = RD(1, 3) RD(3, 2) = RD(2, 3) RD(3, 3) = (D(1, 1) + D(2, 2)  2 * D(1, 2)  2 * D(3, 3)) * Sin(mpa) ^ 2 * Cos(mpa) ^ 2 + D(3, 3) * (Cos(mpa) ^ 4 + Sin(mpa) ^ 4) End Sub '********************************************************************************** Sub StiffnessElement() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 Erase ke For z1 = 1 To 2 For z2 = 1 To 2 t = gp(z1) s = gp(z2) 'Derivatives of shape functions sf1t = (1 + s) / 4 sf1s = (1 + t) / 4 sf2t = (1  s) / 4 sf2s = (1  t) / 4 sf3t = (1  s) / 4 sf3s = (1  t) / 4 sf4t = (1 + s) / 4 sf4s = (1 + t) / 4 'Calculate Jacobian J11 = x1 * sf1t + x2 * sf2t + x3 * sf3t + x4 * sf4t J12 = y1 * sf1t + y2 * sf2t + y3 * sf3t + y4 * sf4t J21 = x1 * sf1s + x2 * sf2s + x3 * sf3s + x4 * sf4s J22 = y1 * sf1s + y2 * sf2s + y3 * sf3s + y4 * sf4s detJ = J11 * J22  J12 * J21 'linear calculate strain gradient matrix BL(1, 1) = (J12 * sf1s + J22 * sf1t) / detJ BL(3, 1) = (J11 * sf1s  J21 * sf1t) / detJ BL(2, 2) = BL(3, 1) BL(3, 2) = BL(1, 1) BL(1, 3) = (J12 * sf2s + J22 * sf2t) / detJ BL(3, 3) = (J11 * sf2s  J21 * sf2t) / detJ BL(2, 4) = BL(3, 3) BL(3, 4) = BL(1, 3) BL(1, 5) = (J12 * sf3s + J22 * sf3t) / detJ BL(3, 5) = (J11 * sf3s  J21 * sf3t) / detJ BL(2, 6) = BL(3, 5) BL(3, 6) = BL(1, 5) BL(1, 7) = (J12 * sf4s + J22 * sf4t) / detJ BL(3, 7) = (J11 * sf4s  J21 * sf4t) / detJ BL(2, 8) = BL(3, 7) BL(3, 8) = BL(1, 7) For m1 = 1 To 8 For m2 = m1 To 8 For m3 = 1 To 3 For m4 = 1 To 3 ke(m1, m2) = ke(m1, m2) + BL(m3, m1) * RD(m3, m4) * BL(m4, m2) * te * detJ Next Next Next Next Next Next End Sub '********************************************************************************** Sub StiffnessConstraint() '*********Apply constraints to Global Stiffness matrix by Elimination approach 'Create dummy stiffness matrix for elimination of constrained DOF ReDim Kdummy(2 * (nw + nc + 1) * (nt + 1), 2 * (nw + nc + 1) * (nt + 1)) For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1) Kdummy(q1, q2) = KG(q1, q2) Next Next 'Identify dof for applying geometric constraints For q3 = 1 To (nc + nw + 1) 'Horizontal axes constraints Yconstraints = 2 * ((nt + 1) * (q3  1) + 1) 'Vertical axes constraints Xconstraints = 2 * ((nt + 1) * (q3))  1 '*********Eliminate constrained dof (delete corresponding rows and columns) XCons = Xconstraints  2 * (q3  1) YCons = Yconstraints  2 * (q3  1) For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) 'Column Eliminations If q2 < YCons And q2 < XCons Then Kdummy(q1, q2) = Kdummy(q1, q2) ElseIf q2 >= YCons And q2 < XCons Then If q2 + 1 = XCons Then Kdummy(q1, q2) = Kdummy(q1, q2 + 2) Else Kdummy(q1, q2) = Kdummy(q1, q2 + 1) End If ElseIf q2 >= XCons Then Kdummy(q1, q2) = Kdummy(q1, q2 + 2) End If Next Next For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) 'Row Eliminations If q1 < YCons And q1 < XCons Then Kdummy(q1, q2) = Kdummy(q1, q2) ElseIf q1 >= YCons And q1 < XCons Then If q1 + 1 = XCons Then Kdummy(q1, q2) = Kdummy(q1 + 2, q2) Else Kdummy(q1, q2) = Kdummy(q1 + 1, q2) End If ElseIf q1 >= XCons Then Kdummy(q1, q2) = Kdummy(q1 + 2, q2) End If Next Next Next '********Calculate constrained global stiffness matrix ReDim Kconstrained(((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))), (2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To ((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q2 = 1 To ((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) Kconstrained(q1, q2) = Kdummy(q1, q2) Next Next 'Worksheets("Kconstrained").Activate 'Worksheets("Kconstrained").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'Range(Cells(5 + q1, 5 + q2), Cells(5 + q1, 5 + q2)) = Kconstrained(q1, q2) 'Next 'Next End Sub '********************************************************************************** Sub MassConstraint() '*********Apply constraints to Global Lumped Mass matrix by Elimination approach 'Create dummy Lumped Mass matrix for elimination of constrained DOF ReDim Mdummy(2 * (nw + nc + 1) * (nt + 1), 2 * (nw + nc + 1) * (nt + 1)) For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1) Mdummy(q1, q2) = MG(q1, q2) Next Next 'Identify dof for applying geometric constraints For q3 = 1 To (nc + nw + 1) 'Horizontal axes constraints Yconstraints = 2 * ((nt + 1) * (q3  1) + 1) 'Vertical axes constraints Xconstraints = 2 * ((nt + 1) * (q3))  1 '*********Eliminate constrained dof (delete corresponding rows and columns) XCons = Xconstraints  2 * (q3  1) YCons = Yconstraints  2 * (q3  1) For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) 'Column Eliminations If q2 < YCons And q2 < XCons Then Mdummy(q1, q2) = Mdummy(q1, q2) ElseIf q2 >= YCons And q2 < XCons Then If q2 + 1 = XCons Then Mdummy(q1, q2) = Mdummy(q1, q2 + 2) Else Mdummy(q1, q2) = Mdummy(q1, q2 + 1) End If ElseIf q2 >= XCons Then Mdummy(q1, q2) = Mdummy(q1, q2 + 2) End If Next Next For q1 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) For q2 = 1 To 2 * (nw + nc + 1) * (nt + 1)  (2 * q3) 'Row Eliminations If q1 < YCons And q1 < XCons Then Mdummy(q1, q2) = Mdummy(q1, q2) ElseIf q1 >= YCons And q1 < XCons Then If q1 + 1 = XCons Then Mdummy(q1, q2) = Mdummy(q1 + 2, q2) Else Mdummy(q1, q2) = Mdummy(q1 + 1, q2) End If ElseIf q1 >= XCons Then Mdummy(q1, q2) = Mdummy(q1 + 2, q2) End If Next Next Next '********Calculate constrained global stiffness matrix ReDim Mconstrained(((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))), (2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To ((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q2 = 1 To ((2 * (nw + nc + 1) * (nt + 1))  (2 * (nc + nw + 1))) Mconstrained(q1, q2) = Mdummy(q1, q2) Next Next 'Worksheets("Mconstrained").Activate 'Worksheets("Mconstrained").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'Range(Cells(5 + q1, 5 + q2), Cells(5 + q1, 5 + q2)) = Mconstrained(q1, q2) 'Next 'Next End Sub '********************************************************************************** Sub PolarPointer() 'uses dummies q1dummy,q2dummy 'call with q1,q2 loop p(1) = 2 * ((q1  1) * (nt + 1) + q2)  1 p(2) = 2 * ((q1  1) * (nt + 1) + q2) p(3) = 2 * ((q1  1) * (nt + 1) + q2 + 1)  1 p(4) = 2 * ((q1  1) * (nt + 1) + q2 + 1) p(5) = 2 * ((q1) * (nt + 1) + q2)  1 p(6) = 2 * ((q1) * (nt + 1) + q2) p(7) = 2 * ((q1) * (nt + 1) + q2 + 1)  1 p(8) = 2 * ((q1) * (nt + 1) + q2 + 1) End Sub '********************************************************************************** Sub PolarWebMat() Erase D matcon = Etweb  Erweb * vrtweb ^ 2 D(1, 1) = Erweb * Etweb / matcon D(1, 2) = Erweb * Etweb * vrtweb / matcon D(2, 1) = D(1, 2) D(2, 2) = Etweb ^ 2 / matcon D(3, 3) = 2 * Erweb ' * Etweb / ((Erweb + Etweb) * (1 + vrtweb)) End Sub '********************************************************************************** Sub MassSystemLumped() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 'clear system arrays ReDim MG(2 * (nc + nw + 1) * (nt + 1), 2 * (nc + nw + 1) * (nt + 1)) '*****Dummy loop for filling 0 values (only for checking) For q1 = 1 To 2 * (nc + nw + 1) * (nt + 1) For q2 = 1 To 2 * (nc + nw + 1) * (nt + 1) MG(q1, q2) = 0 Next Next '**************************************** 'form core region lumped mass matrix For q1 = 1 To nc For q2 = 1 To nt x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) Call MassCoreElementLumped Call PolarPointer For m1 = 1 To 8 For m2 = m1 To 8 MG(p(m1), p(m2)) = MG(p(m1), p(m2)) + Mele(m1, m2) Next Next Next Next 'form web region lumped mass matrix For q1 = nc + 1 To nc + nw For q2 = 1 To nt x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) Call MassWebElementLumped Call PolarPointer 'Assemble Global Lumped Mass matrix (Diagnol matrix) For m1 = 1 To 8 For m2 = m1 To 8 MG(p(m1), p(m2)) = MG(p(m1), p(m2)) + Mele(m1, m2) Next Next Next Next 'Worksheets("MG").Activate 'Worksheets("MG").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'Range(Cells(5 + q1, 5 + q1), Cells(5 + q1, 5 + q1)) = MG(q1, q1) 'Next End Sub '********************************************************************************** Sub MassCoreElementLumped() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 Erase Mele For z1 = 1 To 2 For z2 = 1 To 2 t = gp(z1) s = gp(z2) 'Derivatives of shape functions sf1t = (1 + s) / 4 sf1s = (1 + t) / 4 sf2t = (1  s) / 4 sf2s = (1  t) / 4 sf3t = (1  s) / 4 sf3s = (1  t) / 4 sf4t = (1 + s) / 4 sf4s = (1 + t) / 4 'Calculate Jacobian J11 = x1 * sf1t + x2 * sf2t + x3 * sf3t + x4 * sf4t J12 = y1 * sf1t + y2 * sf2t + y3 * sf3t + y4 * sf4t J21 = x1 * sf1s + x2 * sf2s + x3 * sf3s + x4 * sf4s J22 = y1 * sf1s + y2 * sf2s + y3 * sf3s + y4 * sf4s detJ = J11 * J22  J12 * J21 'Calculate elemntal lumped mass matrix for core region For m1 = 1 To 8 Mele(m1, m1) = (te * detJ * rhocore) / (4 * 386) Next Next Next End Sub '********************************************************************************** Sub MassWebElementLumped() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 Erase Mele For z1 = 1 To 2 For z2 = 1 To 2 t = gp(z1) s = gp(z2) 'Derivatives of shape functions sf1t = (1 + s) / 4 sf1s = (1 + t) / 4 sf2t = (1  s) / 4 sf2s = (1  t) / 4 sf3t = (1  s) / 4 sf3s = (1  t) / 4 sf4t = (1 + s) / 4 sf4s = (1 + t) / 4 'Calculate Jacobian J11 = x1 * sf1t + x2 * sf2t + x3 * sf3t + x4 * sf4t J12 = y1 * sf1t + y2 * sf2t + y3 * sf3t + y4 * sf4t J21 = x1 * sf1s + x2 * sf2s + x3 * sf3s + x4 * sf4s J22 = y1 * sf1s + y2 * sf2s + y3 * sf3s + y4 * sf4s detJ = J11 * J22  J12 * J21 'Calculate elemntal lumped mass matrix for core region For m1 = 1 To 8 Mele(m1, m1) = (te * detJ * rhoweb) / (4 * 386) Next Next Next End Sub '********************************************************************************** Sub MassEffective() '****************initial impact on tip node meff = 2 * ((nt + 1) * (nc + nw + 1)) MG(meff, meff) = MG(meff, meff) + (0.00071456) 'Worksheets("MG").Activate 'Worksheets("MG").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'Range(Cells(5 + q1, 5 + q2), Cells(5 + q1, 5 + q2)) = MG(q1, q2) 'Next 'Next End Sub '********************************************************************************** Sub MassLumpedInverse() ReDim MGINV((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)), (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MGINV(q1, q1) = 1 / Mconstrained(q1, q1) Next 'Worksheets("Minv").Activate 'Worksheets("Minv").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'Range(Cells(5 + q1, 5 + q1), Cells(5 + q1, 5 + q1)) = MGINV(q1, q1) 'Next End Sub '********************************************************************************** Sub MinvF() 'Calculate [MGinv] * [F] ReDim MinF((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MinF(q1) = MinF(q1) + (MGINV(q1, q1) * F0(q1)) Next End Sub '********************************************************************************** Sub MinvMU00() 'Calculate [MinM] * [U00] ReDim MinM((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)), (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MinM(q1, q1) = MinM(q1, q1) + (MGINV(q1, q1) * MG(q1, q1)) Next ReDim MinMU00((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MinMU00(q1) = MinMU00(q1) + (MinM(q1, q1) * U00(q1)) Next End Sub '********************************************************************************** Sub TIBracketU0() 'Calculate 2*[M] for bracketed term in deformation eq ReDim MTI2((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)), (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) MTI2(q1, q1) = MTI2(q1, q1) + (Mconstrained(q1, q1) * 2) Next 'Multiply (2*[M]*d0) ReDim M2d0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) M2d0(q1) = M2d0(q1) + (MTI2(q1, q1) * U0(q1)) Next 'Calculate delT^2*[K] for bracketed term in deformation eq Call SystemStiffnessUpdater Call StiffnessConstraint ReDim KT2((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)), (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) KT2(q1, q2) = KT2(q1, q2) + (Kconstrained(q1, q2) * delT ^ 2) Next Next 'Multiply (delT^2*[K]*d0) ReDim KT2d0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) KT2d0(q1) = KT2d0(q1) + (KT2(q1, q2) * U0(q2)) Next Next 'Calculate bracketed term {2[M](delT^2)[K]}d0 in deformation eq ReDim BrTI((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) BrTI(q1) = BrTI(q1) + (M2d0(q1)  KT2d0(q1)) Next 'Multiply bracketed term by M Inverse {[Minv]*(2[M](delT^2)[K]})*d0)} ReDim BrMinv((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) BrMinv(q1) = BrMinv(q1) + (MGINV(q1, q1) * BrTI(q1)) Next End Sub '********************************************************************************** Sub TimestepAbaqus() ' Calculate smallest element dimension in the mesh (Differenece between Y coordinate of tip node and its adjacent node) t1 = ((nt + 1) * (nc + nw)) + (nt + 1) t2 = ((nt + 1) * (nc + nw  1)) + (nt + 1) Lmin = Abs(YC(t1)  YC(t2)) 'Calculate dilation wave speed in terms of Lame's constant and material density lamda = (Etweb * vrtweb) / ((1 + vrtweb) * (1  (2 * vrtweb))) mu = Etweb / (2 * (1 + vrtweb)) Cd = Sqr((lamda + (2 * mu)) / (rhoweb / 386)) 'Calculate time step delT = ((Lmin / Cd)) * (0.707) * MultiFactor Worksheets("FE Nodes").Activate Worksheets("FE Nodes").Cells.ClearContents 'For q1 = 1 To NN 'Range(Cells(q1 + 5, 5), Cells(q1 + 5, 5)) = XC(q1) 'Range(Cells(q1 + 5, 6), Cells(q1 + 5, 6)) = YC(q1) 'Next End Sub '********************************************************************************** Sub TransientDI() cons = delT nmax = 25000 ReDim U00((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), U0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim U1((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), Un((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim V00((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), V0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim V1((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), Vn((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim A00((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), A0((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim A1((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))), An((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim F0(2 * (nc + nw + 1) * (nt + 1)) ReDim Ag(nmax) '*****Erase worksheet data for previous runs Worksheets("Plot").Activate Worksheets("Plot").Cells.ClearContents 'Worksheets("Deformation").Activate 'Worksheets("Deformation").Cells.ClearContents 'Worksheets("Deceleration").Activate 'Worksheets("Deceleration").Cells.ClearContents 'Input initial vectors For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) U0(q1) = 0 'Displacement vector at t=0 V0(q1) = 0 'Velocity vector at t=0 A0(q1) = 0 'Acceleration vector at t=0 F0(q1) = 0 'External force vector at t=0 Next 'Input initial velocity at t=0 invel = (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) V0(invel) = V0(invel)  (32.31) 'Calculate displacment for time step before t=0 For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) U00(q1) = U00(q1)  (delT * V0(q1)) Next Call TIBracketU0 'Call MinvF 'Call MinvMU00 '******Calculate for nmax time steps For n1 = 1 To nmax ReDim Un((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) ReDim An((2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1))) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) '***Un(q1) = Un(q1) + ((MinF(q1) * delT ^ 2) + (((MinMU0(q1) * 2)  (HU0(q1) * delT ^ 2))  (MinMU00(q1)))) Un(q1) = Un(q1) + ((BrMinv(q1))  (U00(q1))) Next Call Udumstress 'Calculate deceleration from deformations For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'An(q1) = An(q1) + ((HUn(q1))) An(q1) = An(q1) + ((Un(q1)  (2 * U0(q1)) + (U00(q1))) / delT ^ 2) Next tipnode = (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) Ag(n1) = An(tipnode) / 386 '****Print results on worksheets for all nodes 'Worksheets("Plot").Activate 'Range(Cells(5 + n1, 3), Cells(5 + n1, 3)) = delT 'Range(Cells(5 + n1, 4), Cells(5 + n1, 4)) = cons 'Range(Cells(5 + n1, 5), Cells(5 + n1, 5)) = Un(tipnode 'Range(Cells(4 + n1, 7), Cells(4 + n1, 7)) = Ag(n1) 'Range(Cells(4 + n1, 9), Cells(4 + n1, 9)) = Rhos '****Reassign values to vectors for carrying out calculations for next time step For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) U00(q1) = U0(q1) U0(q1) = Un(q1) Next Call CartesianMatUpdater Call TIBracketU0 cons = cons + delT Next 'Calculate hardness from maximum deceleration Rhos = 0 For q1 = 1 To nmax If Ag(q1) > Rhos Then Rhos = Ag(q1) End If Next Rhos = Rhos / 2.6 Worksheets("INPUT").Activate Range(Cells(30, 2), Cells(30, 2)) = Rhos End Sub '********************************************************************************** Sub Udumstress() ReDim Udummy(2 * (nw + nc + 1) * (nt + 1)) '2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1) ReDim Ucopy(2 * (nw + nc + 1) * (nt + 1)) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) Ucopy(q1) = Un(q1) Next 'Identify dof for applying geometric constraints For q3 = 1 To (nc + nw + 1) 'Horizontal axes constraints Yconstraints = 2 * ((nt + 1) * (q3  1) + 1) 'Vertical axes constraints Xconstraints = 2 * ((nt + 1) * (q3))  1 '*********Eliminate constrained dof (delete corresponding rows and columns) XCons = Xconstraints '+ 2 * (q3  1) YCons = Yconstraints '+ 2 * (q3  1) For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) + (2 * q3) 'Row addition: Add zero deformation to constrained DOF If q1 < YCons And q1 < XCons Then Udummy(q1) = Ucopy(q1) ElseIf q1 = YCons Then Udummy(q1) = 0 ElseIf q1 > YCons And q1 < XCons Then Udummy(q1) = Ucopy(q1  1) ElseIf q1 = XCons Then Udummy(q1) = 0 ElseIf q1 > XCons Then Udummy(q1) = Ucopy(q1  2) End If Next '*****Update Ucopy vector For q2 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) + (2 * q3) Ucopy(q2) = Udummy(q2) Next Next 'Worksheets("Dummy").Activate 'Worksheets("Dummy").Cells.ClearContents 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1)) 'Range(Cells(5 + q1, 4), Cells(5 + q1, 4)) = Udummy(q1) 'Next 'For q1 = 1 To (2 * (nc + nw + 1) * (nt + 1))  (2 * (nc + nw + 1)) 'Range(Cells(5 + q1, 2), Cells(5 + q1, 2)) = Un(q1) 'Next End Sub '********************************************************************************** Sub CartesianMatUpdater() ReDim PRStress((nc + nw) * (nt)) For q1 = nc + 1 To nc + nw For q2 = 1 To nt NE = (q1  1) * nt + q2 x1 = XC((q1  1) * (nt + 1) + q2) x2 = XC((q1  1) * (nt + 1) + q2 + 1) x3 = XC((q1) * (nt + 1) + q2) x4 = XC((q1) * (nt + 1) + q2 + 1) y1 = YC((q1  1) * (nt + 1) + q2) y2 = YC((q1  1) * (nt + 1) + q2 + 1) y3 = YC((q1) * (nt + 1) + q2) y4 = YC((q1) * (nt + 1) + q2 + 1) mpa = (THETAC(q2 + 1) + THETAC(q2)) / 2 'mean polar angle Erweb = K2 * (K1  PRTotal(NE)) Call PolarWebMat Call CartesianMatCon Call PolarPointer Call CartesianStress '************add stress increment to initial stress PRStress(NE) = PRStress(NE) + (cstress(1) * Cos(mpa) ^ 2 + cstress(2) * Sin(mpa) ^ 2 + 2 * cstress(3) * Cos(mpa) * Sin(mpa)) / 4 Next Next End Sub '********************************************************************************** Sub CartesianStress() 'dummy variables 'x1,x2,x3,x4,y1,y2,y3,y4 Erase cstress 'erase cartesian stress components For z1 = 1 To 2 For z2 = 1 To 2 t = gp(z1) s = gp(z2) 'derivatives of shape functions sf1t = (1 + s) / 4 sf1s = (1 + t) / 4 sf2t = (1  s) / 4 sf2s = (1  t) / 4 sf3t = (1  s) / 4 sf3s = (1  t) / 4 sf4t = (1 + s) / 4 sf4s = (1 + t) / 4 'calculate jakobien J11 = x1 * sf1t + x2 * sf2t + x3 * sf3t + x4 * sf4t J12 = y1 * sf1t + y2 * sf2t + y3 * sf3t + y4 * sf4t J21 = x1 * sf1s + x2 * sf2s + x3 * sf3s + x4 * sf4s J22 = y1 * sf1s + y2 * sf2s + y3 * sf3s + y4 * sf4s detJ = J11 * J22  J12 * J21 'linear calculate strain gradient matrix BL(1, 1) = (J12 * sf1s + J22 * sf1t) / detJ BL(3, 1) = (J11 * sf1s  J21 * sf1t) / detJ BL(2, 2) = BL(3, 1) BL(3, 2) = BL(1, 1) BL(1, 3) = (J12 * sf2s + J22 * sf2t) / detJ BL(3, 3) = (J11 * sf2s  J21 * sf2t) / detJ BL(2, 4) = BL(3, 3) BL(3, 4) = BL(1, 3) BL(1, 5) = (J12 * sf3s + J22 * sf3t) / detJ BL(3, 5) = (J11 * sf3s  J21 * sf3t) / detJ BL(2, 6) = BL(3, 5) BL(3, 6) = BL(1, 5) BL(1, 7) = (J12 * sf4s + J22 * sf4t) / detJ BL(3, 7) = (J11 * sf4s  J21 * sf4t) / detJ BL(2, 8) = BL(3, 7) BL(3, 8) = BL(1, 7) For m1 = 1 To 3 For m2 = 1 To 3 For m3 = 1 To 8 cstress(m1) = cstress(m1) + RD(m1, m2) * BL(m2, m3) * Udummy(p(m3)) Next Next Next Next Next End Sub '********************************************************************************** ADVISER‟S APPROVAL: Dr. J. K. Good VITA Dipesh Dilip Mistry Candidate for the Degree of Master of Science Thesis: CONNECTING WINDING CODES TO QUALITY CONTROL DEVICES. Major Field: Mechanical and Aerospace Engineering Biographical: Personal Data: Born in Vasai, India 28th September 1985, son of Mr. Dilip D. Mistry and Mrs. Darshana D. Mistry. Education: Received Diploma in Mechanical Engineering from Bhausaheb Vartak Polytechnic, Vasai in June 2004. Received Bachelor of Science degree from University of Pune, Pune, in May 2007. Completed the requirements for the Master of Science majoring in Mechanical and Aerospace Engineering at Oklahoma State University, Stillwater, Oklahoma in July 2010. Experience: Worked as a teaching assistant for Advanced Experimental Methods in Design (Jan 2008 to May 2010) and Finite Element Methods (Aug 2009 to Dec 2009) at OSU. Worked as a Graduate Research Assistant at Web handling Research Center (WHRC), OSU (June 2008 to July 2010). ADVISER‟S APPROVAL: Dr. J. K. Good Name: Dipesh Dilip Mistry Date of Degree: July, 2010 Institution: Oklahoma State University Location: Stillwater, Oklahoma Title of Study: CONNECTING WINDING CODES TO QUALITY CONTROL DEVICES. Pages in Study: 62 Candidate for the Degree of Master of Science Major Field: Mechanical and Aerospace Engineering Scope and Method of Study: The goal of this study was to develop an extension to the winding models which can convert the outputs of pressures and stresses to units that can be measured using existing quality control devices, the Rho meter in this case. The Rho meter makes dynamic measurement of the hardness of a wound roll. The peak deceleration of the Rho meter striker is measured and converted into arbitrary “Rho” units, which is a measure of hardness of the wound roll. Initially, experiments were conducted to study the variation of hardness in a wound roll as a function of winding tension and the roll outer radius. Based on the experimental results, a dynamic impact model was developed to simulate the Rho meter. The output obtained from the model was compared to the experimental results. Findings and Conclusions: The model made fair predictions of the roll hardness. The material properties vary in time due to the pressures generated by the impact. The explicit method was used in the dynamic analysis since it facilitates the updating of material properties after each solution in time step. Due to this, large computation times are associated with this model. Different algorithms were used to calculate the time step required for the explicit solution and the algorithm employed in this model shows good convergence. The accuracy of the model output increases with an increase in mesh density. ADVISER‟S APPROVAL: Dr. J. K. Good 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


