UFDC Home  myUFDC Home  Help 



Full Text  
FINITE ELEMENT ANALYSIS FOR INCIPIENT FLOW OF BULK SOLID IN A DIAMONDBACK HOPPER By OSAMA SULEIMAN SAADA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2005 Copyright 2005 by Osama Suleiman Saada This document is dedicated to: my father Suleiman Saada, my mother Nema Aude and the rest of my family members for without their support this research would not have been completed. ACKNOWLEDGMENTS I would like to thank the following, because without their help this research would not have been possible: my committee chairman, Dr. Nicolaie Cristescu, for his continuous help; Dr. Kerry Johanson, who was very kind to advise and guide me throughout my graduate research; Dr. Ray Bucklin for his support, not only in scientific decisions but also in academic and personal ones; the rest of my committee members Dr. Frank Townsend and Dr. Bhavani Sankar. I would also like to extend my deepest gratitude to another member of my committee, Dr. Olesya Zhupanska, who gave me detailed help on the modeling part of the research. I would also like to thank the Particle Engineering Center for funding this project. Special thanks go to the PERC faculty and staff for their assistance. Last but not least I would like to extend special thanks to my father, Suleiman Saada, and my mother, Nema Aude, for their patience and continuous moral and financial support for the many years of my college career. TABLE OF CONTENTS ACKNOW LEDGM ENTS ........................................ iv LIST OF TABLES ...... ......... ..... .... ........................... vii LIST OF FIGU RE S ........................................... ................................... viii ABSTRACT.................. .................. xii CHAPTER 1 INTRODUCTION ................... ............................ ......... .. .......... 1 2 BACKGROUND INFORM ATION ........................................................................4 B ulk Solid C classification .........................................................................4 Flow Patterns in Bins/Hoppers ..........................................4 The Diam ondback H opper ............................................... ........ 10 Discharge Aids......................................................... .... ......... ...............1 Janssen or Slice M odel A analysis ........................................ ................. 12 Silo Vertical Section......................................... .........12 Converging H opper .......................................................... .......... ....... ........14 Diam ondback HopperTM ....................................................... 15 Direct Shear TestersJenike Test ................................ ............ ............18 Jenike W all Friction Test....................................................... 21 Schulze Test....................................................................22 Indirect Shear TesterTriaxial Test.............. ................................ ....... 23 The Diamondback HopperTM Measurements .................................. .....28 3 FINITE ELEM EN T ................................. .....................................33 Background in Finite Elem ent M odeling ..................................... ......... ......33 FEM using ABAQUS on Diamondback................................. ..........................36 Explicit Time Integration ................. ...............................37 Geom etry, M eshing, and Loading ....................... ........................... ...............43 Plasticity M odels: General D iscussion..............................................................49 v 4 CONSTITUTIVE MATERIAL LAWS ..........................................55 Capped DruckerPrager M odel ................................. .................... .... ........... 58 3D General Elastic/Viscoplastic Model................... ............................................59 Numerical Integration of the Elastic/Viscoplastic Equation...............................63 5 EXPERIMENTAL RESULTS AND DISCUSSION...............................................67 Param eter D term inationShear Tests ...................................................... 67 Diamondback hopper measurementsTekscan Pads.......................73 Conclusions and Discussion of the Test results and Testing Procedure...................77 6 FEM RESULTS AND DISCUSSIONS ....................................................79 Capped DruckerPrager M odel ................................. .................... .... ........... 79 Drucker Prager M odel Verification .............................. ............... 81 Drucker Prager Model Hopper FEM Results...........................82 Conclusions And Discussion of the Predictive Capabilities of the Capped D ruckerPrager M odel .................................................. ............... 91 Viscoplastic Model Parameter Determination.................................................93 Time Effects ...................................... ....... ........... .98 M odel V alidation ..................................... ......................... ..... .............. 101 Viscoplastic Model Hopper FEM Results ...................... .....103 Conclusions and Discussion of the Predictive Capabilities of the 3D Elastic/Viscoplastic Model ..................... ......... ..................108 7 CONCLUSIONS ................ .... ............. ............... ...10 APPENDIX A NUMERICAL INTEGRATION SCHEME FOR TRANSIENT CREEP................113 B TE ST IN G PR O C E D U R E .................................................................................... 118 Sam ple Preparation ................... ............................. ......... .... ..... ................. 118 M embrane Correction............................................. 120 Piston Friction Error .................. .............................................. .. ........ 121 C TEK SC A N D A TA ................................................................................. 122 D DRUCKERPRAGER FEM DATA............................ ......... 124 E VISCOPLASTIC M ODEL FLOW CHART...........................................................129 L IST O F R E F E R E N C E S ............................................................................................ 130 BIOGRAPH ICAL SKETCH .............. ................................................................. .....135 LIST OF TABLES Table page 2.1 Classification of powders ......................... ...... ........4 2.2 Comparison of mass flow and funnel flow of particulate materials .......................6 2.3 Pad specifications ................ ........ .. ............. ........ .... .. .... .. 30 3.1 Numerical simulation research projects ..................... ............ 34 5.1 Summary of Schulze test results ................................. ............... 68 5.2 B elt velocity vs. flow rate ............................................... ............... 73 6.1 Capped DruckerPrager parameters determined from shear tests............................81 6.2 Coefficients of the yield function.................................................................. 95 6.3 Coefficients of the viscoplastic potential .............................................................. 97 LIST OF FIGURES Figure page 2.1 Types of flow obstruction ................................................. ...............5 2.2 Patterns offlow ........................................... .........6 2.3 Types of hoppers ................ ..... ... ....... .. .... ........ .............7 2.4 Problem s in silos ............................................... .......... 8 2.5 Trajectories of m ajor principle stress ...................... ........... ................................9 2.6 Schematic of the Diamondback HopperTM showing position of pad .....................10 2.7 Forces acting on differential slice of fill in bin ....................................................13 2.8 Forces acting in converging hopper .......................................................................14 2.9 Stresses acting on differential slice element in diamondback hopperTM...............16 2.10 Flow wall stresses at an axial position Z=3.8 cm below the hopper transition........ 17 2.10 Yield loci of different bulk solids .................................. ...........18 2.13 Yield loci .......................................................20 2.14 Jenike w all friction test ............................................................... 21 2.15 Schulze rotational tester .............................................................................22 2.16 E xperim ental setup..................................................23 2.17 Triaxial chamber on an axial loading device ................. ................. ...........24 2.18 Filling procedure ....................................................... ...............25 2.19 Stressstrain curves .............................................. ..... ...27 2.20 Layers of discs dilating as they are sheared .................................. .....28 2.21 TekScanTM Pads ................. ..............................29 2 .22 T ek Scan P ads ....................................................... 30 2.23 Location of TekscanTM pads in the diamondback hopper.................... ...........31 3.1 Flowchart showing steps followed to complete an analysis in ABAQUS............37 3.2 Summary of the explicit dynamics algorithm ...................................................42 3.3 Geometry and meshing of hopper wall ..............................43 3.4 Filling procedure for powder...................... ........ ...............44 3.5 Geometry and meshing of bulk solid .......................................... 46 3.6 Frictional behavior ......................................... ...... .... .....48 4.1 Tresca and MohrCoulomb yield surfaces .........................................55 4.2 Von M ises and DruckerPrager yield surfaces.............................. ........... ....56 4.3 Closed yield surface ...................... .......... ........ ........ 57 4.4 The linear DruckerPrager cap model ........................................... 58 4.5 Domains of compressibility and dilatancy ..............................................62 5.1 Output data from Schulze test on Silica at 8Kg ........................................... 67 5.2 Schulze test results for to determine the linear DruckerPrager surface..................68 5.3 Jenike wall friction results to determine boundary conditions...............................69 5.4 Hydrostatic triaxial testing on Silica Powder (5 confining pressures)...................70 5.5 Axial deformation of 4 deviatoric triaxial tests on Silica Powder and a rate of 0.1 of mm/min ............... ............... ................. 71 5.6 Volumetric deformation of 4 deviatoric triaxial tests on Silica Powder and a rate of 0. 1mm/min ............... .................................... ..71 5.7 Axial deformation of 4 deviatoric tests on Silica Powder at 13.8KPa and 3 4 5 K p a ...................... .. .. ......... .. .. ....................... 7 2 5.8 D ischarge m echanism and belt.................................... .................. 73 5.9 Pad location and output example ................................. ............... 74 5.10 Static Wall stress data at a distance 2.8cm below hopper transition........................75 5.11 Wall stress data at a distance 2.8cm below hopper transition at various speeds and static conditions ................... ........................................................................ ..........75 5.12 Wall stress data at a distance 5.6cm below hopper transition at various speeds and static conditions ................... ........................................................................ ..........76 5.13 Fourier series analysis on the wall pressure measurements. .............................77 6.1 The Linear DruckerPrager Cap model........................................ 79 6.2 FEM simulation of a triaxial test using the DruckerPrager Model.........................81 6.3 Axial stressstrain curves for experiment and FEM using the DruckerPrager m odel ................... ......... .......................... .......................... 82 6.4 FEM calculated contact stresses at various filling steps .......................................83 6.5 FEM calculated contact stresses at various discharge steps .....................................84 6.6 FEM vs. TekScan area of interest ............................... ............... 84 6.7 FEM calculated static wall stresses vs. measured wall stresses at seven locations ........................................................86 6.8 FEM calculated contact stresses vs. measured wall stresses at seven locations during flow ...................................... ................................. ......... 87 6.9 FEM calculated contact stresses vs. measured wall stresses at two locations during flow at two different times........................ ........ ............................ 88 6.10 FEM calculated contact stresses vs. measured wall stresses at two locations during flow angles of friction values.................................................89 6.11 FEM calculated contact stresses vs. measured wall stresses at two locations during flow with various cohesion values...............................................................90 6.12 The irreversible volumetric stress work (data points) and function HH(solid line) ................... ............ .. ....................94 6.13 The irreversible volumetric stress work (data points) and function HD(solid line) ................... ............ .. ....................94 6.14 The viscoplastic potential derivatives. ....................................... 97 6.15 History dependence of the stabilization boundary. ...............................................98 6.16 The irreversible volumetric stress work at two rates 0. 1mm/min and 1mm/min.....99 6.17 Deviatoric creep test at a confining pressure of 34.5 Kpa .....................................100 6.18 Variation in time of axial creep rate at a confining pressure of 34.5 Kpa..............100 6.19 Variation in time of volumetric creep rate at a confining pressure of 34.5Kpa ................... ... ............................ ......... 101 6.20 Theoretically predicted stressstrain curves (solid lines) vs. experimentally determined curves at 34.5 Kpa confining pressure .................................... 103 6.21 FEM calculated static contact stresses vs. measured wall stresses at seven locations below the hopper transition............... .................................... 104 6.23 FEM calculated flow contact stresses vs. measured wall stresses at two locations below the hopper transition at two wall friction values ......................................107 6.24 FEM calculated flow contact stresses vs. measured wall stresses at two locations below the hopper transition using adaptive meshing ..........................................107 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy FINITE ELEMENT ANALYSIS FOR INCIPIENT FLOW OF BULK SOLID IN A DIAMONDBACK HOPPER By Osama Suleiman Saada December 2005 Chair: Nicolaie Cristescu Major Department: Mechanical and Aerospace Engineering Storing of bulk materials is essential in a large number of industries. Powders are often stored in containers such as a bin, hopper or silo and discharged through an opening at the bottom of the container under the influence of gravity. This research tackles the frequent problems encountered in handling bulk solids such as flow obstructions and discontinuous flow resulting in doing and piping in hoppers. This problem is of interest to a variety of industries such as chemical processing, food, detergents, ceramics and pharmaceuticals. An objective of this work is to use finite element analysis to produce results that could easily be incorporated into the industry and be used by the design engineers. The study includes both a numerical approach, with an appropriate material response model, and an experimental procedure implemented on the bulk material of interest. Results from the analysis are used to relate the slopes of the channels and the size of the outlets. Results from the FEM analysis are verified using measured stresses from a DiamondbackTM hopper. CHAPTER 1 INTRODUCTION Almost every industry handles powders/bulk solids, either as raw materials or as final products. Examples include chemicals and chemical processing, foods, detergents, ceramics, minerals, and pharmaceuticals. Powders are often stored in a container such as a bin, hopper or silo and removed through an opening in the bottom of that container under the influence of gravity. The reliability of the process involved depends on the flowability of these powders. Most powders are cohesive, that is, have the tendency to agglomerate or stick together over time. For the material to flow out of a storage facility, bridging, arching or doing must be prevented. For a stable arch to form, the bulk solid must gain enough strength to support itself within the constraints of the container. The strength is a function of the degree of the compaction of the material and stress is a function of spatial position in a piece of process equipment. Thus knowing the stress allows us to compute the strength of the bulk solid. Stresses are also effective in producing yield or failure of the bulk solid. Blockage occurs when the strength exceeds the stresses needed to fail the material. Since the advent of more powerful computers, numerical methods have become very useful in research on flow of bulk solids. Numerical methods are very economical and lend themselves to comprehensive parametric studies. This study presents a finite element approach to solve for displacements, velocities and stresses of a cohesive powder. The domain is discretized into small elements and displacements and loads are approximated. The commercially available ABAQUS software is used. Both the Drucker Prager yield criterion and a 3D viscoplastic model are used to describe the material response. To verify the numerical results, wall pressure measurements inside a diamondback hopperTM are taken. Measurement of wall loads in hoppers has been difficult due to the expense of constructing multiple loads cell measurement units in a hopper. Recent advances in measurement techniques have allowed significant improvements in the measurement of multiple point normal stresses in hoppers. This dissertation presents measurements of wall stresses using pressure sensitive pads made by TekScanTM Reliable and complete data on the deformation, failure and flow behavior will allow a fundamental understanding of powder flow initiation and flow. An objective of this work is to produce results that could easily be incorporated into the industry and be used by the design engineers. Both the finite element approach, with an appropriate material response model, and the technique for measurements of wall stresses on a diamondback hopperTM using pressure sensitive pads are tools towards achieving that objective. Results from the analysis are used to relate the slopes of the channels and the size of the outlets necessary to maintain the flow of a solid of given flowability on walls of given frictional properties. Using the correct material response models, the right boundary conditions combined with the correct testing procedure under the correct loading conditions, continuum models with improved predictive capabilities can be developed. The understanding of powder behavior as well as the amount of information on material properties is limited. Generally, the material properties are determined based on shear test measurements, yet in order to predict the bulk mechanical response during storage and transport, experimental data corresponding to a variety of deformation, loading conditions, and shear rates are needed. To fully understand the material behavior of bulk solids, experiments are required where all possible stress or strain cases are induced in the material specimen. However, this would be an impossible task and simplified tests are instead utilized. In these tests only some stress or strain state components vary independently. In powder mechanics and technology the most commonly used testing devices are shear testers, such as the Jenike shear cell. This type of device has some limitations such as the fact that the location of the shear plane in the sample is determined by the shape of the tester and that the bulk solid is assumed to obey a rigidhardening/softening plastic behavior of MohrCoulomb type. A combination of direct shear testers (the Jenike tester and the Schulze tester), and the indirect shear tester (such as the triaxial tester) with finite element techniques will improve the analysis techniques of powder flow. The specific objectives of this research are: * Present measurements of wall stresses on a diamondback hopperTM using pressure sensitive pads made by TekScanTM * Present a finite element approach to solve for displacements, velocities and stresses of a bulk material based on a viscoplastic model * Verify the method by comparing FEM results with measured stresses on a diamondback hopperTM * Experimentally determine the material properties needed to use in the constitutive relationships * Use the commercially available FEM software ABAQUS with built in models for material response such as the DruckerPrager model * Write a user defined subroutine in FORTRAN to be used with ABAQUS to describe the viscoplastic behavior. CHAPTER 2 BACKGROUND INFORMATION Bulk Solid Classification A particulate system containing particles of about 100[lm and smaller is called a powder. If the majority of the particles in the system are larger in size, it will be called a granular material. A classification of particulate solids based on particle size is given in table 2.1(Brown and Richards, 1970 and Nedderman, 1992). Granular materials are generally free flowing. To initiate flow in granular material it is sufficient to overcome their friction resistance. Most fine powders are prone to agglomerate and stick together, causing significant storage and handling problems Table 2.1: Classification of powders Classification Particle size range(pm) Ultrafine powder <1 Superfine powder 110 Granular powder 10100 Granular solid 1003000 Broken solid >3000 Source: Brown and Richards, 1970 and Nedderman, 1992 Flow Patterns in Bins/Hoppers Before proceeding with a description of the flow patterns in storage facilities the definition of terms that are used throughout this study is necessary. "Hoppers" have inclined sidewalls while "Bins" is the term used to describe a combination of a hopper and a vertical section. Silos is also a term used to describe a tall vertical container where the height to diameter ratio is bigger than 1.5. In European literature this term is used to describe a vertical section on top of a hopper section. These containers vary in their usage in industrial applications. In agricultural, mining, cement and refractory applications they are used as "primary" storage facilities and have big capacities (up to 1000 tons). In other industrial applications these containers are used only for short time storage, the bulk material being subsequently transferred to other containers or mixers. Thus, they are generally much smaller than the "primary hoppers." In bins, most powders can experience obstruction to flow as shown in figure 2.1. Doming Doing Piping Figure 2.1: Types of flow obstruction For an arch to form, the solid needs to have developed enough strength to support the weight of the obstruction. Some of the methods used to reinitiate flow are vibration of the bin and manual movement of the powder. These methods have proven to be costly, inefficient and dangerous. There are two major kinds of flow patterns in a container: (1) mass flow, and (2) funnel flow. In mass flow, the entire material flows throughout the whole vessel during discharge. Figure 2.2a illustrates such flow. In funnel flow, there is a stagnant layer of material at the wall of the vessel. This stagnant region extends all the way up to the top of the container. Flow from a funnel flow bin occurs by first emptying the center flow channel and then, provided the material is sufficiently free flowing, material sloughs off the top surface. If the material is cohesive a stable rathole can form. Flow patterns are illustrated in Figure 2.2b and 2.2c. a) Mass flow b) Semimass flow c) Funnel flow Figure 2.2: Patterns of flow For most applications, the ideal flow situation is mass flow and it occurs when the hopper walls are steep and smooth enough and when there are no abrupt transitions between the bin section and the hopper section. When the walls are rough and the slope angle is too large the flow regime tends to be a funnel flow. Table 2.2 lists the advantages and disadvantages of both mass and funnel flow. Table 2.2: Comparison of mass flow and funnel flow of particulate materials Mass flow Funnel flow Characteristics No stagnant zones Stagnant zone formation Uses full crosssection of vessel Flow occurs within a portion of vessel crosssection Firstin, firstout flow Firstin, lastout flow Advantages Often minimizes segregation, agglomeration Small stresses on vessel walls during flow due to of materials during discharge the 'buffer effect' of stagnant zones. Very low particle velocities close to vessel walls; reduced particle attrition and wall wear Disadvantages Large stresses on vessel walls during flow Promotes segregation and agglomeration during flow Attrition of particles and erosion and wear Discharge rate less predictable as flow region boundary Of vessel wall surface due to high particle can alter with time Velocities Small storage volume/vessel height ratio Source: Brown and Richards, 1970 Storage containers can have axisymmetric or planar geometry as seen in figure 2.3. The shape, number and position of the discharge outlet can vary depending upon the silo geometry, bulk solid properties and process requirements. The geometry of the container (C) CYLINDRICAL FLATBOTTOMED SLOT OPENING r Massflow hoppers E)WEDGE PYRAMID Funnelflow hoppers (D) CYLINDRICAL FLATBOTTOMED CIRCULAR OPENIF Figure 2.3: Types of hoppers and properties of the powder determine the flow pattern during filling and discharge. In many industrial applications the same containers are used for storing different types of bulk solid. These materials have different mechanical and physical properties such as particle size, particle distribution, particle shape, and bulk density and frictional properties. This results in varying flow regimes for different materials. A hopper designed to provide mass flow for one material may not provide mass flow for a different material. Figure 2.4 shows the relationships among bulk solid material properties, the silo filling process, the flow patterns during discharge, the wall pressures, the stress induced in the structure and the conditions that might cause failure of the structure (Rotter 1998). (B) SQUARE OPENING (C) CHISEL A) CONICAL HOPPER F AMID. i i RE OPENING BQ (B) CONICAL The flow pattern is strongly affected by the distribution of densities and orientation of particles which are determined by the mechanical characteristics of the solid and the Aspect of silo behavior Loss of Function Solid properties ARCHING/RATHOL Filling Methods '' SEGREGATION Flow pattern Pressure on silo walls Stresses in silo structure Failure conditions for the structure  COLLAPSE Figure 2.4: Problems in silos way it is filled. These characteristics also determine whether arching or ratholing occurs at the outlet. The flow pattern during discharge also determines if segregation occurs and influences the pressure on the walls. Hence in order to predict the nonuniform and unsymmetrical wall pressure it is essential that the flow pattern be predicted. This influences directly the wall stress conditions that could induce failure and collapse of the structure. In regards to bulk solids handling in storage facilities, 3 phases can be distinguished: (1)filling, (2)storage and (3)discharge. In the filling stage, with the outlet closed, the material is supplied slowly at the top of the container. As more and more material is poured in, stresses, due to the weight of layers, start developing throughout the container. In the storage stage no material is flowing in or out of the container. For cohesive (fine) powders this may lead to an increase of interparticle bonding. Due to its own weight the stored material compacts and its strength increases with time. Also, the slow escape of the air in pores results in stronger surface contact and an increase in frictional forces. In the third and final stage the material is discharged through an opening at the bottom of the container. Upon opening the outlet, a dilation zone neighboring the outlet is formed. This zone will spread to the upper layers of the material and cause it to flow. In hoppers with steep and smooth sides, this zone of dilated material will cover the whole crosssectional area of the container resulting in massflow, while for rough walls the dilation zone is confined only to the center, resulting in funnel flow. The state of stress throughout the container, in the discharge stage, drastically differs from the stress state in the filling/storage phases. The corresponding stress states are called active and passive, respectively. In the filling and storage stages the direction of the major principal stress c0I is vertical with slight curvature close to the walls. In the discharge stage the principal stress is horizontal with slight curvature at the walls (Figure 2.5). Active Passive a) Filling/storage b) discharge Figure 2.5: Trajectories of major principle stress The Diamondback HopperTM The Diamondback HopperTM consists of a roundtooval hopper with diverging end walls. This hopper section is positioned above an ovaltoround hopper that necks down to a circular outlet. Pressure sensitive pads (from TekScanTM) were used to measure the normal wall pressure (Johanson, 2001). These pads consist of two plastic sheets with conductive pressure sensitive material printed in rows on one sheet of plastic and columns on the other sheet of plastic. When these sheets contact each other they form conductive junctions where contact resistance varies with the normal stress applied. The effective area of each junction is 1 cm2 and there are over 2000 independent normal stress measurements possible that can be recorded at cycle times up to eight per second. The voltage drops across these junctions were sequentially measured and scaled to give real engineering force units using the data acquisition software. The TekScanTM pad was glued to the inside surface of the roundtooval hopper with one edge of the pad along the center of the flat plate section (see Figure 2.6). Pad 61 cm Pad locations 2064 sensors Figure 2.6: Schematic of the Diamondback HopperTM showing position of pad Discharge Aids Poor design or incorrect use of a hopper leads to the use of devices that reduce flow problems. The use of the wrong device, however, may have the reverse effects and create more problems than it solves. These devices are used if it proves that there are constraints in the overall system that prevent the unaided gravity flow of the material. These aids developed from practices such as beating the hopper with a blunt instrument and poking the material with some sort of rod. The three major types of aids are (1) pneumatic relying on the application of air to the product; (2) vibrationalrelying on mechanical vibration of the hopper or the material; (3) mechanicalphysically extracting the product from the hopper. In the pneumatic aids, aeration devices are used to introduce air at the time that the material is discharged so as to "fluidize" the material in the region of the outlet opening and to reduce the friction between the material and the hopper wall and the second is to introduce a "trickle of flow" of air during the whole period that the product is stored to prevent the gain of strength in the material. Air is also introduced into inflatable pads that act mechanically against the stored material. When using vibrational methods, devices could be used to vibrate the hopper or bin walls or to vibrate the material directly. Vibration should not be applied when the outlet is closed, as this could result in the strengthening of any arch formation. In mechanical methods powered dislodgers such as vertical or horizontal stirrers are used to manually move the material . Janssen or Slice Model Analysis Silo Vertical Section Janssen (1895) developed a method of predicting vertical silo wall pressures over a century ago that is still widely used in the industry although most of the assumptions he used in the derivations have been shown to be incorrect. The simplest slice for the analysis is the planar finite boundary slice that has sides that are perpendicular to the walls of the bin. The forces acting on an element slice of stored powder, of thickness dz, at a depth h from the top surface in a deep bin with an overall height h, crosssectional area A and circumference C are shown in figure 2.7. Let the vertical stress at the upper surface, depth h, be ov and at the lower surface, depth h+dz, be ov+dov. ow and Tw are the normal and tangential (shear) stresses at the wall due to friction between the walls and the bulk material. The weight of the powder in the slice is Adzpg where g is the acceleration due to gravity and p is the material bulk density which is assumed to remain constant over the entire depth of the powder. Equating the forces in the vertical direction we get UA + Adzpg A(a, + do,) + r,Cdz = 0 (2.1) It is assumed that the ratio of horizontal to vertical pressure is constant everywhere in the element: U = K (2.2) 0v For fully mobilized friction we have /=tan W(2.3) OyA h rCdz ___ dz (oC+dov)A Adzpg Figure 2.7: Forces acting on differential slice of fill in bin Expanding equation 2.1 and simplifying we get the first order differential equation: Separating variables and integrating with the boundary condition ov=0 at h=0 (at the level of the free surface), and ov=ov at h=h gives gA kPC h ua =KC e )A (2.4) and pgA k h] a = Ka 1 e (2.5) For large depth of fill the maximum values of vertical and normal stresses are oa (max) =pgA (2.6) puKC and a, (max) =gA (2.7) tiC while near the top surface:  (top) pgh (2.8) The ratio A/C varies for different silo geometries (e.g., rectangular, circular). As mentioned above recent studies have challenged the accuracy of Janssen's assumptions of constant K values, angle of wall friction and bulk density. The value of K has a considerable influence on the stress distribution in the material and on the wall. Researchers have differing views of the K value. Jenike (1973), for example, assumed K=0.4 for most granular materials while others claim it is a function of internal angle of friction. Converging Hopper The above analysis is applicable to deep silos and flat bottom bins ignoring any end effects. For the hopper section a linear hydrostatic pressure gradient is proposed 9, = pg[ o +h] (2.9) P9 where avo is the vertical pressure at the top of the hopper calculated from equation 2.5 as shown in figure 2.8. \ zcyw I Figure 2.8: Forces acting in converging hopper Walker (1966) assumed a constant value of K equal to: K = tan (2.10) (tan p, + tan a) While Jenike (1961) assumed a constant value of: sin (2a cos) (2.11) [sin (, + 2a) + sin 9] Diamondback HopperT A force balance can be done on a small differential slice of material in the Diamondback HopperTM as shown in Figure 2.9 (Johanson and Bucklin, 2004). The forces acting on the material slice element are due to the bulk material slice weight, wall frictional forces on the flat plate section and the round end walls, vertical stresses, and normal stress conditions at the bin wall. The definition of Kvalue and the columbic friction condition are used to relate the vertical pressure acting on the material to the stresses normal and tangent to the bin wall (see Equations 2.12 through 2.15). The resulting differential equation is found in Equation 2.16. The crosssectional area (A) and hopper dimensions (L) and (D) are functions of the axial coordinate as indicated in Equations 2.16 through 2.21. L = KL .v (2.12) rL = KL .a,tan(w) (2.13) a = K, (a) c, (2.14) z, = Ke (a) a, tan(O) (2.15) dz Lr I^It'^ a Figure 2.9: Stresses acting on differential slice element in diamondback hopperTM 1 dA 2K, (LD)) dA + 2.K (L D)(tan(O) + tan(D))+ di A dz A =Yg 4.D (2.16) A K, (a) (tan() + tan(O(a))) da 0 A=D +(LD)D (2.17) d= 2F [ D+L Ltan(OD)2Dtan(OL) (2.18) L L 2.ztan(OL) (2.19) D= D 2 z tan(OD) (2.20) K, (a) = (K1m, Klmax) (sin(a))" +K1ma (2.21) Where values of parameters Kimin, Kimax, and n were chosen to minimize the deviation between the observed wall stress profiles and the calculated wall stresses over the entire axial hopper depth. The hopper wall slope angle changes around the perimeter. The slope angle along the flat plate section equals the hopper angle (OD) while the slope angle along the hopper end wall equals the hopper angle (OL). The variation in the hopper angle, around the perimeter, is given by Equation 2.23, where. qa is the hopper section angle as defined in Figure 2.9. 0(a) = a tan[( sin(a)). tan(OD) + sin(a) tan(OL)] (2.22) Figure 2.10 shows Janssen computed wall stresses compared with the actual TekscanTM measurements for the section of the hopper shown (Johanson, 2003). 7 4 Hopper Section Ange, (deg) Computed Janssen Stress A Fully Developed Flow Figure 2.10: Flow wall stresses at an axial position Z=3.8 cm below the hopper transition 1 :A 0 10 20 30 40 50 60 70 80 90 Hopper Section Angle, a (deg) Computed Janssen Stress A Fully Developed Flow Figure 2. 10: Flow wall stresses at an axial position Z=3.8 cm below the hopper transition The Janssen slice model can provide a first approximation to the loads in diamondback hoppers. However, there is significant variation from the simple slice model approach used here. More complex constitutive models will be required to increase the agreement between experimental and theoretical approaches. Direct Shear TestersJenike Test For the past 30 years the assessment of flow properties of bulk solids has been done through shear testing (Jenike, 1961). The bulk solid is assumed to obey a rigid hardening/softening plastic behavior of MohrCoulomb type (Figure 2.10) for coarse cohesionless, particles (granular particles), the yield locus is approximated by: cT=atan4 (2.23) while for fine powders, the yield locus is c=atan4 + c (2.24) a) Yield locus of a b) Yield locus of a c) Yield locus of a free cohesive Coulomb cohesive non solid Coulomb solid Figure 2.10: Yield loci of different bulk solids As shown in figure 2.11, the cell consists of a base, an upper ring and a cover. The cover has a fixture that can hold a normal force N. The motor causes the loading stem to push against the shear ring and displace it horizontally. The shear force is measured and is divided by the crosssectional area of the ring to provide the shear stress(c). SHEAR FORCE X HEAR FORCE POWDER SAMPLE Figure 2.11: Jenike shear test The test begins when the sample of bulk solid is poured into the shear cell and is preconsolidated manually by placing a special top on the cell, applying a normal load and twisting the top to consolidate the material. Then a normal stress on is applied and the specimen is sheared. The shear force continuously increases with time until constant shear is obtained. Corresponding to the change in shear is a change in bulk density. After a period of time, the bulk density reaches a constant value for a value of on. At that point, the deformation is known as steady state deformation. The sample is first presheared under a constant on and T. The shearing is stopped when steadystate deformation is reached. Then the sample is sheared under the normal stress on that is smaller than the shear stress. To insure the same initial bulk density the sample is then presheared under the same normal stress on and sheared at a lower on than the first test. This process is repeated several times. The yield locus, usually called the effective yield locus is used to determine the flow parameters. As seen from figure 2.12 it is a straight line passing through the origin and tangential to the Mohrs circle at steady state deformation (the largest circle). This line has the inclination 4e or the effective angle of internal friction (angle with the o axis). At lower normal stresses the yield loci is curved instead of straight. Linearized Yield locus Yield locus .,* . a e Figure 2.12: Yield loci Therefore a linearized yield locus is used to define the internal friction angle (4i). Using a different on for preshearing the sample will result in different consolidation and bulk density values. This will produce several yield loci and several Mohr circles as seen in Figure 2.13 FC Fc 21 Y Figure 2.13: Yield loci Jenike Wall Friction Test The Jenike tester (1961) can also be utilized to determine the wall friction angle. The test is modified so that the base of the tester is replaced by a flat plate of the wall material that the hopper is made of (Figure 2.14). APPLIED WEIGHT '.NORMAL FORCE) SHEAR FORCE POWDER SAMPLE Figure. 2.14: Jenike wall friction test The procedure is similar to the description above in the sample preparation stage. A maximum normal force is applied at the top cover and is decreased in a series of steps while the maximum shearing force is being measured. As explained above the wall yield locus can be plotted with c vs o and the coulomb equation is: = arctan (u) (2.25) The coefficient of friction is expressed as the angle of wall friction given by: T = Pu( (2.26) BRACKET SHEAR FORCE Schulze Test The Jenike test is limited to a small displacement of about 6mm and is suitable for fine particles only. Rotational shear testers such as the annular Schulze (1994) tester (figure 2.15) have unlimited strains. The powder is contained in the cell and loaded from F, a F2 Figure 2.15: Schulze rotational tester the top with a normal force N through the lid. After calibrating and filling the cell with powder a predetermined consolidation weight was added to a hanger, which is connected to the top of the cell. The machine was then turned on and allowed to run until it reaches a steady state. The normal weight was adjusted to about 70% of the steady state value and the sample was sheared. During the test the shear cell rotates slowly in the direction of o while the cover is prevented from rotation by two tie rods. This causes the powder to shear. The forces Fi and F2 and the normal forces are recorded. This procedure was repeated several times with decreasing weight. Between each test the sample was pre consolidated to the steady state value. The data reduction procedure is same as the Jenike test and gives the yield loci and the internal angle of friction. Indirect Shear TesterTriaxial Test To investigate the stressstrain behavior of powders, a triaxial testing apparatus is used. Its capabilities are wider than those of uniaxial testers, and it can be described as an axially symmetric device having 2 degrees of freedom. A schematic diagram and picture of the triaxial apparatus are shown in Figure 2.16 and 2.17. Water Rubber membrane Volume Change Device Porous stone * Figure 2.16: Experimental setup In an ideal triaxial test, all three major principal stresses would be independently controlled. However the independent control would lead to mechanical difficulties that limit the conduction of such tests to special applications. The commonly used triaxial test refers to the axisymmetric compression test. This test has been used since the beginning of the 19th century in testing the strength of soils. The soil sample is precompacted, saturated with water, and compressed both hydrostatically and deviatorically and the water discharged from the sample is used to measure the volume change and hence the density change. In recent years, the test has been used to test other bulk solids: grains and fine powders. The testing procedure was modified so as to allow testing of dry powders at low stresses. The behavior in this regime is essential for the applications of interest. Figure 2.17: Triaxial chamber on an axial loading device A latex membrane is stretched out into a cylindrical shape using a 2piece cylindrical mold and vacuum pressure. The powder is poured into the membrane dispersed and compressed using a predetermined stress. Figure 2.18 shows the steps followed in the filling procedure. The value of AH was determined using Janssen's equation. To insure repeatable initial conditions, both the mass and volume of the material have to be measured. The membrane is sealed by rubber Orings to a pedestal at the bottom, and to a cap at the top. Vacuum grease is used the membrane and the pistons so as to reduce the chances of leaks. Any test consists of two phases: (1) a hydrostatic phase followed by (2) a deviatoric phase. In the hydrostatic test the assembly is placed in a chamber that is filled with water. The water is pressurized to the desired confining pressure. This hydrostatic pressure is increased up to a desired value ca using a pressurizing device. The sample can support itself due to the filling procedure and thus the mould is removed. Under the applied hydrostatic pressure the powder, compacts and the specimen takes the form of a cylinder. The change in volume of the specimen is recorded using a data acquisition system that measures the volume of the water displaced. The initial volume of the sample is determined by multiplying the height by the arithmetic mean of the diameters at 3 locations in the vicinity of the middle of the sample. Figure 2.18: Filling procedure In the deviatoric phase, the assembly is placed on an axial loading device. An additional axial stress is applied by means of a piston passing through a frictionless I~ bushing at the top of the chamber. The test is carried out under strain control conditions in which a predetermined rate of axial deformation is imposed and the axial load required to maintain the rate of deformation is measured via a load cell. As the sample deforms, water is displaced from the chamber. The quantities measured during the deviatoric phase are: the confining pressure 03 (which provides an allaround pressure on the lateral surface of the sample), the axial force applied to the piston, the change in length of the sample and the change in the volume of the sample which is determined by the change in the volume of water existing in the chamber. Thus, in the test, the change in volume of the sample we can be monitored continuously as it is sheared. These tests will produce the following stressstrain curves: the axial strain E, versus the deviatoric stress a, Ca (o, is the axial stress) and the volumetric strain Es versus the deviatoric stress a7 1 3, respectively. Here and throughout the text compressive stresses and strains are considered to be positive. The axial strain is defined as: ei= where 41 is the initial length of the sample 10 and / is the current length ; the volumetric strain is E, = Vo Vo and V being the 0V initial and current volume, respectively. Any change in the volumetric strain reflects change of the relative positions of the powder particles. Failure corresponds to an observed instability in the specimen when one part slides with respect to the other along a plane inclined at about 450 with respect to the vertical axis of the specimen. Failure shows up in the stress strain curve as a maximum deviatoric pressure. 300 4* S250 p 150 elv(test 3) 100 * 50 Powder Test Date u3(kPa) pl(g/cm^3) Poly#1 3 12/10 4056 08594 55000 35000 15000 5000 25000 45000 65000 85000 105000 Srain(E6) Figure 2.19 Stressstrain curves Figure 2.19 shows the stressstrain curves obtained in a typical test (Saada, 1999). This is a test carried out on a powder at a confining pressure of 40.56kPa. The axial strain curve shows that the powder undergoes elastic deformation up to about 50% of its strength, followed by work hardening and failure. The stressvolumetric strain curve shows two distinct regimes of behavior: first, the powder compacts (i.e. Sv increases) very slightly, then starts to expands (i.e. S; decreases). The volume increase, which is observed on the last portion of the curve (o, c3) Sv is called dilatancy. The powder expands because of the interlocking of the particles: deformation can proceed only if some particles are able to ride up or rotate over other particles. This phenomena can simplistically be visualized with the aid of Figure 2.20, which shows two layers of discs, one on top of the other. If a shear stress is applied to the upper layer, then each disc in this layer has to rise (increasing the volume occupied by the particulate material) for the sample to undergo any shear deformation. Figure 2.20: Layers of discs dilating as they are sheared The Diamondback HopperTM Measurements The core of the TekscanTM pressure measurement system (1987) consists of an extremely thin 0.004 in (0.1 mm), flexible tactile force sensor. Sensors come in both grid based and single load cell configurations, and are available in a wide range of shapes, sizes and spatial resolutions (sensor spacing). These sensors are capable of measuring pressures ranging from 015 kPa to 0175 MPa. Each application requires an optimal match between the dimensional characteristics of the objects) to be tested and the spatial resolution and pressure range provided by Tekscan's sensor technology. Sensing locations within a matrix can be as small as .0009 square inches (.140 mm2); therefore, a one square centimeter area can contain an array of 170 of these locations. Teksan's Virtual System Architecture (VSA) allows the user to integrate several sensors into a uniform whole. The standard sensor consists of two thin, flexible polyester sheets which have electrically conductive electrodes deposited in varying patterns as seen in Figure 2.21 (Tekscan Documentation, 1987). In a simplified example below, the inside surface of one sheet forms a row pattern while the inner surface of the other employs a column pattern. The spacing between the rows and columns varies according to sensor application and can be as small as 0.5 mm. Figure 2.21: TekScanTM Pads Before assembly, a patented, thin semiconductive coating (ink) is applied as an intermediate layer between the electrical contacts (rows and columns). This ink, unique to Tekscan sensors, provides the electrical resistance change at each of the intersecting points. When the two polyester sheets are placed on top of each other, a grid pattern is formed, creating a sensing location at each intersection. By measuring the changes in current flow at each intersection point, the applied force distribution pattern can be measured and displayed on the computer screen. With the TekscanTM system, force measurements can be made either statically or dynamically and the information can be seen as graphically informative 2D or 3D displays. In use, the sensor is installed between two mating surfaces. Tekscan's matrixbased systems provide an array of force sensitive cells that enable measurement of the pressure distribution between the two surfaces (figure 2.22). The 2D and 3D displays show the location and magnitude of the forces exerted on the surface of the sensor at each sensing location. Force and pressure changes can be observed, measured, recorded, and analyzed throughout the test, providing a powerful engineering tool. Overall Width (W) I Matrix Width (MW) * Tab Length (A) U ____ Matrix Height .IH, SRow Width (rw) Column Width (cw) Su Row Spacing (rs) Column Spacing (cs) Sensell Exploded View Figure 2.22: TekScan Pads Table 2.3: Pad specifications General Dimensions Sensing Region Dimensions Summary Overall Overall Tab Matrix Matrix Length Width Length Width Height Columns Rows # of Sensels L W A MW MH CW CS Qty RW RS Qty Sensels Density (mm) (mm) (mm) (mm) (mm) (mm) (mm) (mm) (mm) (per cm2) 622 530 130 488 427 6.35 10.2 48 6.35 10. 2 42 2016 0.97 Several pleats or folds were made in the upper portion of the pad to make it conform to the upper cylinder wall (figure 2.23). A thin retaining ring was placed at the transition between the hopper and cylinder. This held the pad close to the hopper wall Overall Length (L) Pad locations 2064 sensors Figure 2.23: Location of TekscanTM pads in the diamondback hopper and helped maintain the pleats in the cylinder. Contact paper was placed on the hopper and cylinder surface to protect the pad and produce a consistent wall friction angle (of about 17 degrees) in the bin. The pad lead was fed through the hopper wall well above the transition and connected to the data acquisition system. Spot calibration checks were made on groups of four load sensors. This arrangement produced a load measurement system capable of measuring normal wall loads to within +10 %. The bin level was maintained by a choke fed standpipe located at the bin centerline and at an axial position about 60 cm above the hopper transition. The material used was fine 50 micron silica. During normal operation, the feed system maintains a repose angle at the bottom of the standpipe and produces a relatively constant level as material discharges from the hopper. Flow from the hopper was controlled by means of a belt feeder below the Diamondback HopperTM. There was a gate valve between the Diamondback HopperTM and the belt. Initially, this gate was closed and the lower hopper was charged with a small quantity of material. The amount of material in the hopper initially was enough to fill only the lower ovaltoround section of the bin. The gate was then opened, allowing material to fill the void region between the gate and the belt. This was done to avoid surges onto the belt during the initial filling routine. The remaining hopper was then filled slowly, using the existing conveying system until the diamondback test hopperTM and a surge bin above this test hopper were filled. Material level was maintained in the surge hopper by recycling any material leaving the belt into the surge hopper feeding the Diamondback HopperTM bin. Once the bin was full, material was allowed to stand at rest for about 10 minutes as material deaerated. The speed of the belt was preset and flow was initiated. The wall load data acquisition system recorded changes in wall stress as flow was initiated. The initial stress condition shows a large peak pressure concentrating at the top point of the triangular flat plate in the hopper. This peak is focused at the tip of the plate. During flow these loads decrease and de localize producing the highest peak stresses at along the edge of the plate. CHAPTER 3 FINITE ELEMENT Background in Finite Element Modeling Great effort has been made in recent decades to understand flow phenomena of bulk solids. Since the 1970's, a large number of research teams have worked on the application of finite element analysis to hopper problems. Models and programs were hampered due to limited capacity of computers and the high cost of equipment. Nowadays, it is possible to use computers whose capacity and speed is continually increasing and a large number of programs exist that manipulate, analyze and present the results. In principle there are no restrictions with regards to the hopper geometry and to the bulk material in FEM analysis. However studies have to account for complex geometry of the container, complex material behavior and interaction with the wall. Hopper geometry is not always axisymmetric and often has shapes that require threedimensional numerical modeling, as is the case with the diamondback hopper. Bulk solids exhibit complex mechanical behavior such as anisotropy, plasticity, dilatancy and so on. Most of these properties are developed during discharge due to the large strains that are experienced. Some properties, such as anisotropy are developed during the filling process. It is very essential that detailed mathematical models, accompanied with a proper experimental procedure, are used to describe these features. The appropriate simulation has to be used to model the interaction between the solid and the wall. For the case of smooth wall, a constant friction coefficient can be satisfactory while for rough walls it might be required to use more sophisticated modeling. In a book edited by C.J. Brown and J. Nielsen (1998) several FEM research on silo flow are presented. This book is a compilation of reports brought about by a research project on silos called the CASilo Project. It was funded by the European commission and completed in 1997. The first chairman of the group, G. Rombach (1998), assembled data for comparison of existing programs. Table 3.1 lists several projects regarding granular solid behavior simulation. It includes both finite element and discrete element simulations. Table 3.1: Numerical simulation research projects Name of Static/ 2D/3D Cohesive/ Author/user program Dyna Non mic Cohesive 1. Aubry (Ecole Central Paris) GEFDYN S/D 2D/3D C/N 2. Eibl (Karlsruhe U, Germany) SILO S/D 2D/3D C/N 3. Eibl (Karlsruhe U, Germany) S/D 2D/3D C/N 4. Eibl (Karlsruhe U, Germany) ABAQUS S/D 2D/3D C/N 5. Klisinski (Lulea U, Sweden) SILO S/D 2D N 6. Klisinski (Lulea U, Sweden) BULKFEM S/D 2D N 7. Martinez (INSA, Rennes, France) AMG/PLAXIS S/D 2D N 8. Martinez (INSA, Rennes, France) MODACSIL S 2D/3D C/N 9. Schwedes (Braunschweig U, Germany) ABAQUS S/D 2D/3D C/N 10. Schwedes (Braunschweig U, Germany) HAUFWERK S 2D N 11. Tuzun (University of Surrey, UK) HOPFLO S/D 2D/3D N 12. Ooi, Carter (Edinburgh U, Scotland) AFENA S 2D C/N 13. Ooi (Edinburgh U, Scotland) ABAQUS S/D 2D/3D C/N 14. Rong (Edinburgh U, Scotland) DEM.F S/D 2D N 15. Thompson (Edinburgh U, Scotland) SIMULEX S/D 2D N 16. Cundall (INSA, Rennes, France) TRUBALL,PFC S/D 2D/3D N Also listed in the table are: the name of the program used, whether it analyzed static or dynamic material, the geometry used (2D or 3D) and whether it analyzed cohesive or noncohesive material. The material models used include the wellknown elasticplastic models and yield criteria such as MohrCoulomb, Tresca and DruckerPrager. Other law include: Lade (1977), Kolymbas (1998), Boyce (1980), Wilde (1979) and critical state. These laws have been modified for cohesive bulk solids and dynamic calculations. Aside from this project there has not been extensive work published on finite element modeling flow of cohesive material in hoppers. Haussler and Eibl (1984) developed a model using 18 governing equations: 3 for dynamic equilibrium, 6 constitutive relations and 9 kinematic relations to describe flow behavior. They used FEM, with triangular elements, to determine the spatial response of sand flow from a hopper. They concluded that the velocity in the bin section was constant, while in the hopper section it was maximum at the center line. Schmidt and Wu (1989) used a similar FEM procedure except that they used a modified version of the viscosity coefficient and used rectangular elements. Runesson and Nilsson (1986) used the same dynamic equations and modeled the flow of granular material as a viscousplastic fluid with a Newtonian part and a plastic part. The DruckerPrager (1952) yield criteria was used to represent purely frictional flow. Link and Elwi (1990) used an elasticperfectly plastic model with wall interface elements to describe incipient flow of cohesionless material. Several steps were used to simulate the filling process in layers while incipient release as opposed to full release was simulated for discharge. Flow was detected by a sudden increase in displacements at the outlet followed by failure of the FEM to converge to a solution. The material was simulated using eightnode isoparametric elements while the interface used sixnode isoperimetric elements. The wall pressure results were close to the results obtained by Jenike while the maximum outlet pressure was lower than the Jenike analysis. The results, however were not verified by experimental procedures. A FEM analysis, with a new secant constitutive relationship, was used, by Meng and Jofriet (1992), to simulate cohessionless granular material (soy bean) flow. The outlet velocities were comparable with preliminary experimental results. The solution converged quickly and was stable throughout the simulation. Another model, developed by Diez and Godoy (1992), used viscoplastic behavior and incompressible flow to describe cohesive materials. This model used the DruckerPrager (1952) yield theory and was applied to conical and wedge shaped hoppers. Results obtained compared well with other published work for the conical hopper for cohessionless materials except at the bottom of the hopper. Kamath and Puri (1999) used the modified Camclay equations in a FEM code to predict incipient flow behavior of wheat flour in a mass flow hopper bin. Incipient flow was characterized by the characterization of the first dynamic arch. This arch represented the transition from the static to the dynamic state during discharge. Incipient flow is assumed to occur when the mesh at the outlet of the hopper displays 7% or more axial strain. The bin was discretized using rectangular 4node elements while the hopper is discritized using quadrilateral elements. A thin wall interface element was the for the powder wall interaction. The FEM results were verified experimentally using a transparent plastic laboratory size mass flow hopper bin. The FEM results were within 95% of the measured values. FEM using ABAQUS on Diamondback The numerical model is based on a consistent continuum mechanics approach. In this project the commercially available software ABAQUS (Inc. 2003) is being used. It is a generalpurpose analysis product that can solve a wide range of linear and nonlinear problems involving the static, dynamic, thermal, and electrical response of components. It has several built in models to describe material response. Such materials include: metals, plastics, concrete, sand and other frictional material. It also allows for a user subroutine to be written that allows any material to be implemented. Some of these plastic models are listed in figure 3.1 (ABAQUS 2003). The figure also shows the steps involved in an analysis. The geometry and meshing and material properties are created in the preprocessor (ABAQUS/CAE), which creates an input file that is submitted to the processor (ABAQUS/Standard) or (ABAQUS/Explicit). The results of the simulation can be viewed as a data file or as visualization in ABAQUS/Viewer. Plasticity Preprocessing User Materials ABAQUS/CAE or other software User Subroutines allows any n Input file: model to be implemented Input file: job.inp Creep Volumetric Swelling Simulation Twolayer viscoplasticity ABAQUS/Standard Or ABAQUS/Explicit Extended DruckerPrager Capped DruckerPrager Output files: cam( job.obd, job.dat MohrCoulomb job.res, job.fil res, oCrushable Foam Postprocessing Strainratedependent Plasticity ABAQUS/CAE or other software Figure 3.1: Flowchart showing steps followed to complete an analysis in ABAQUS Explicit Time Integration Exact solutions of problems where finding stresses and deformations on bodies subjected to loading, requires that both force and moment equilibrium be maintained at all times over any arbitrary volume of the body. The finite element method is based on approximating this equilibrium requirement by replacing it with a weaker requirement, that equilibrium must be maintained in an average sense over a finite number of divisions of the volume of the body. Starting from the force equilibrium equation for the volume: Stds + fdV = 0 (3.1) S V where: V is the volume occupied by a part of the body in the current configuration, S is the surface bounding this volume. t is the force per unit of current f is the body force per unit volume and t = n : (T where c is the Cauchy stress matrix and n is the unit outward normal. From the principles of continuum mechanics the deferential form of the equilibrium equation of motion is derived as: V .07 + f = 0 (3.2) The displacementinterpolation finite element model is developed by writing the equilibrium equations in "weak form". This produces the virtual work equation in the classical form: fu :SDdV = fSvtdS+ fv fdV (3.3) V S V where: 9D is virtual strain rate (virtual rate of deformation) 9v is a virtual velocity field This equation tells us that the rate of work done by the external forces subjected to any virtual velocity field is equal to the rate of work done by the equilibrating stresses on the rate of deformation of the same virtual velocity field. The principle of virtual work is the "weak form" of the equilibrium equations and is used as the basic equilibrium statement for the finite element formulation. ABAQUS can solve problems explicitly or implicitly (using ABAQUS/Standard) For both the explicit and the implicit time integration procedures, the equilibrium is M ii = P I (34) where: P are the external applied forces I are the internal element forces M is the mass matrix ii are the nodal accelerations Both procedures solve for nodal accelerations and use the same element calculations to determine the internal element forces. In the implicit procedure a set of linear equations is solved by a direct solution method to obtain the nodal accelerations. This method however becomes expensive and time consuming compared to using the explicit method. The implicit scheme uses the full Newton's iterative solution method to satisfy dynamic equilibrium at the end of the increment at time t + At and compute displacements at the same time. The time increment, At, is relatively large compared to that used in the explicit method because the implicit scheme is unconditionally stable. For a nonlinear problem each increment typically requires several iterations to obtain a solution within the prescribed tolerances. The iterations continue until several quantitiesforce residual, displacement correction, etc.are within the prescribed tolerances. Since the simulation of flow in a diamondback hopper contains highly discontinuous processes, such as contact and frictional sliding, quadratic convergence may be lost and a large number of iterations may be required. Explicit schemes are often very efficient in solving certain classes of problems that are essentially static and involve complex contact such as forging, rolling, and flow. Flow problems are characterized by large distortions, and contact interaction with the hopper wall. It also has the advantage of requiring much less disk space and memory than implicit simulation. To get accurate solutions for a bulk solid in a hopper a relatively dense mesh with thousands of elements needs to be used. And since the explicit method shows great cost savings over the implicit method it is used in this research. Hence using an explicit scheme where the solution is determined without iterating by explicitly advancing the kinematic state from the previous increment is preferred. A central difference rule, to integrate the equations of motion explicitly through time, using the kinematic conditions at one increment to calculate the kinematic conditions at the next increment, is used. The accelerations at the beginning of the current increment (time t) are calculated as: t) =M (P I) (3.5) Since the explicit procedure always uses a diagonal, or lumped, mass matrix, solving for the accelerations is trivial; there are no simultaneous equations to solve. The acceleration of any node is determined completely by its mass and the net force acting on it, making the nodal calculations very inexpensive. The accelerations are integrated through time using the central difference rule, which calculates the change in velocity assuming that the acceleration is constant. This change in velocity is added to the velocity from the middle of the previous increment to determine the velocities at the middle of the current increment: Atl + A1l)) t+ At ( I t) (3.6) The velocities are integrated through time and added to the displacements at the beginning of the increment to determine the displacements at the end of the increment: u i .\ = u ( \ + At u, st 1 A/ U (t+At) UI(t) + A (t+At) t+ (3.7) Thus, satisfying dynamic equilibrium at the beginning of the increment provides the accelerations. Knowing the accelerations, the velocities and displacements are advanced "explicitly" through time. The term "explicit" refers to the fact that the state at the end of the increment is based solely on the displacements, velocities, and accelerations at the beginning of the increment. This method integrates constant accelerations exactly. For the method to produce accurate results, the time increments must be quite small so that the accelerations are nearly constant during an increment. Since the time increments must be small, analyses typically require many thousands of increments. Fortunately, each increment is inexpensive because there are no simultaneous equations to solve. Most of the computational expense lies in the element calculations to determine the internal forces of the elements acting on the nodes. The element calculations include determining element strains and applying material constitutive relationships (the element stiffness) to determine element stresses and, consequently, internal forces. Here is a summary of the explicit dynamics algorithm: 4 Integrate explicitly through time. (AtI(A) + AtI() uI(i = u 2 t u  u(t+At) u (t) + A (t+At) t+ Element calculations. Compute element strain increments, de, from the strain rate, E Compute stresses, a, from constitutive equations A (tet)= f( o(rcs) Assemble nodal forces, (tAt) Set t + At to t Figure 3.2: Summary of the explicit dynamics algorithm Contact conditions and other extremely discontinuous events are readily formulated in the explicit method and can be enforced on a nodebynode basis without iteration. The nodal accelerations can be adjusted to balance the external and internal forces during contact. The most striking feature of the explicit method is the lack of a global tangent stiffness matrix, which is required with implicit methods. Since the state of the model is advanced explicitly, iterations and tolerances are not required. Nodal calculations. Dynamic equilibrium. id )= M 1 (P I ) Geometry, Meshing, and Loading The advantage of FEM is that a variety of problems such as the hopper geometry, the size and locations of the openings and the interaction between the material and the walls can be investigated. The diamondback hopper geometry is relatively complicated relative to other hopper geometries. Because of this complexity generating the mesh with the appropriate density and correct elements is not a trivial task. The steel wall is simulated as a rigid body as shown in figure 3.3 (from different view points). Since this research focused on the behavior of the bulk material in the hopper, no finite element analysis was conducted on the wall. Because of the complexity of the geometry however, it was required to discritize the wall. It was meshed into 6872 triangular elements (R3D3) as seen in figure 3.3. The mesh had to be of particular density so as to simulate contact analysis correctly. Figure 3.3: Geometry and meshing of hopper wall It was decided to simulate the process of filling the bulk solid in the hopper in several steps. The solid was simulated as a deformable body and divided into 11 different sets. A uniform load is applied to each set to simulate gravity (9.81m/s2) in 11 different steps (as seen in figure 3.4). Discharge was simulated in a 12th step by applying a small displacement to the bottom surface of the powder. By default, all previously defined loads are propagated to the current step. The starting condition for each general step is the ending condition of the previous general step. Thus, the model's response evolves during a sequence of general steps in a simulation. Figure 3.4: Filling procedure for powder It was decided to use second order elements as opposed to first order elements. This is because of the standard firstorder elements are essentially constant strain elements and the solutions they give are generally not accurate and, thus, of little value. The second order elements are capable of representing all possible linear strain fields and much higher solution accuracy per degree of freedom is usually available with the higherorder elements. However, in plasticity problems discontinuities occur in the solution if the finite element solution is to exhibit accuracy, these discontinuities in the gradient field of the solution should be reasonably well modeled. With a fixed mesh that does not use special elements that admit discontinuities in their formulation, this suggests that the firstorder elements are likely to be the most successful, because, for a given number of nodes, they provide the most locations at which some component of the gradient of the solution can be discontinuous (the element edges). Therefore firstorder elements tend to be preferred for plastic cases. Tetrahedral elements are geometrically versatile and are used in many automatic meshing algorithms. It is very convenient to mesh a complex shape with tetrahedra, and the secondorder and modified tetrahedral elements (C3D10, C3D1OM,) in ABAQUS are suitable for general usage. However, tetrahedra are less sensitive to initial element shape. The elements become much less accurate when they are initially distorted. A family of modified 6node triangular and 10node tetrahedral elements is available that provides improved performance over the firstorder tetrahedral elements and that avoids some of the problems that exist for regular secondorder triangular and tetrahedral elements, mainly related to their use in contact problems. The modified tetrahedron elements use a special consistent interpolation scheme for displacement. Displacement degrees of freedom are active at all userdefined nodes. These elements are used in contact simulations because of their excellent contact properties. In ABAQUS/Explicit these modified triangular and tetrahedral elements are the only secondorder elements available. In addition, the regular elements may exhibit volumetricc locking" when incompressibility is approached, such as in problems with a large amount of plastic deformation. The modified elements are more expensive computationally than lowerorder quadrilaterals and hexahedron and sometimes require a more refined mesh for the same level of accuracy. However, in ABAQUS/Explicit they are provided as an attractive alternative to the lowerorder tetrahedron to take advantage of automatic tetrahedral mesh generators. In the diamondback the powder body is meshed into quadratic tetrahedron (C3D10OM) elements and the wall into discrete rigid body. The deformation of the wall is negligible compared to the deformation of the powder. The mesh density and element type are two factors that have a major influence on the accuracy of results and computer time. Hence in this research the appropriate mesh had to be chosen. Figure 3.5 shows 4 different meshes that were used for analysis. As the project progressed, one mesh was used with the various models. Figure 3.5: Geometry and meshing of bulk solid The wall friction coefficient is applied between the wall and the body using Coulomb criteria. The basic concept of the Coulomb friction model is to relate the maximum allowable frictional (shear) stress across an interface to the contact pressure between the contacting bodies. In the basic form of the Coulomb friction model, two contacting surfaces can carry shear stresses up to a certain magnitude across their interface before they start sliding relative to one another; this state is known as sticking. The Coulomb friction model defines this critical shear stress, r,,, at which sliding of the surfaces starts as a fraction of the contact pressure, p between the surfaces ( T, = ,up). The stick/slip calculations determine when a point transitions from sticking to slipping or from slipping to sticking. The fraction, ji, is known as the coefficient of friction. The basic friction model assumes that u is the same in all directions isotropicc friction). For a threedimensional simulation there are two orthogonal components of shear stress, r, and rz, along the interface between the two bodies. These components act in the slip directions for the contact surfaces or contact elements. ABAQUS combines the two shear stress components into an "equivalent shear stress," 7, for the stick/slip calculations, where T = r 2 + r, In addition, ABAQUS combines the two slip velocity components into an equivalent slip rate, Pe 712 +2 The stick/slip calculations define a surface. The friction coefficient is defined as a function of the equivalent slip rate and contact pressure: P = {(eqp ) (3.8) where y,, is the equivalent slip rate, p is the contact. The solid line in Figure 3.6 summarizes the behavior of the Coulomb friction model: there is zero relative motion (slip) of the surfaces when they are sticking (the shear stresses are below IIP)). i (shearstress) slipping sticking Y (slip) Figure 3.6: Frictional behavior Simulating ideal friction behavior can be very difficult; therefore, by default in most cases, ABAQUS uses a penalty friction formulation with an allowable "elastic slip," shown by the dotted line in Figure 125. The "elastic slip" is the small amount of relative motion between the surfaces that occurs when the surfaces should be sticking. ABAQUS automatically chooses the penalty stiffness (the slope of the dotted line) so that this allowable "elastic slip" is a very small fraction of the characteristic element length. The penalty friction formulation works well for most problems, including most metal forming applications. In those problems where the ideal stickslip frictional behavior must be included, the "Lagrange" friction formulation can be used in ABAQUS/Standard and the kinematic friction formulation can be used in ABAQUS/Explicit. The "Lagrange" friction formulation is more expensive in terms of the computer resources used because ABAQUS/Standard uses additional variables for each surface node with frictional contact. In addition, the solution converges more slowly so that additional iterations are usually required. This friction formulation is not discussed in this guide. Kinematic enforcement of the frictional constraints in ABAQUS/Explicit is based on a predictor/corrector algorithm. The force required to maintain a node's position on the opposite surface in the predicted configuration is calculated using the mass associated with the node, the distance the node has slipped, and the time increment. If the shear stress at the node calculated using this force is greater than 7Trit, the surfaces are slipping, and the force corresponding to Tritis applied. In either case the forces result in acceleration corrections tangential to the surface at the slave node and the nodes of the master surface facet that it contacts. Often the friction coefficient at the initiation of slipping from a sticking condition is different from the friction coefficient during established sliding. The former is typically referred to as the static friction coefficient, and the latter is referred to as the kinetic friction coefficient. Plasticity Models: General Discussion The elasticplastic response models in ABAQUS have the same general form. They are written as rateindependent models or as ratedependent. A rateindependent model is one in which the constitutive response does not depend on the rate of deformation as opposed to a ratedependent models where time effect such as creep are considered. A basic assumption of elasticplastic models is that the deformation can be divided into an elastic part and a plastic part. In its most general form this statement is written as: F=Fe.FPl (3.9) where: F is the total deformation gradient, F"1 is the fully recoverable part of the deformation FPl is the plastic deformation The rigid body rotation at the point can be included in the definition of either F"1 or F"1 or can be considered separately before or after either part of the decomposition. This decomposition can be reduced to an additive strain rate decomposition, if the elastic strains are assumed to be infinitesimal: 8 = 8 el + pl (3.10) where: E is the total strain rate s"1 is the elastic strain rate P'1 is the plastic strain rate The strain rate is the rate of deformation: ~Sv~ S= sym (3.11) The above decomposition implies that the elastic response must always be small in problems in which these models are used. In practice this is the case: plasticity models are provided for metals, soils, polymers, crushable foams, and concrete; and in each of these materials it is very unlikely that the elastic strain would ever be larger than a few percent. The elastic part of the response derived from an elastic strain energy density potential, so the stress is defined by: a = (3.12) where: U is the strain energy density potential The stress tensor a is defined as the Cauchy stress tensor. For several of the plasticity models provided in ABAQUS the elasticity is linear, so the strain energy density potential has a very simple form. For the soils model the volumetric elastic strain is proportional to the logarithm of the equivalent pressure stress. The rateindependent plasticity models have a region of purely elastic response. The yield function, f defines the limit to this region of purely elastic response and, for purely elastic response is written so that: f(O*, H,) <0 (3.13) where: 0 is the temperature H, are a set of hardening parameters(the subscript o is introduced simply to indicate that there may be several hardening parameters) In the simplest plasticity model (perfect plasticity) the yield surface acts as a limit surface and there are no hardening parameters at all: no part of the model evolves during the deformation. Complex plasticity models usually include a large number of hardening parameters. Only one is used in the isotropic hardening metal model and in the Camclay model; six are used in the simple kinematic hardening model. When the material is flowing inelastically the inelastic part of the deformation is defined by the flow rule, which we can write as: d c = d (3.14) where ds = is plastic strain increase dA = proportionality factor for the ith system g (, H,,) = plastic potential for the ith system a = Cauchy stress state In an "associated flow" plasticity model the direction of flow is the same as the direction of the outward normal to the yield surface: Df, S= i (3.15) where: c, is a scalar. These flow models are useful for materials in which dislocation motion provides the fundamental mechanisms of plastic flow when there are no sudden changes in the direction of the plastic strain rate at a point. They are generally not accurate for materials in which the inelastic deformation is primarily caused by frictional mechanisms. The metal plasticity models in ABAQUS (except cast iron) and the Camclay soil model use associated flow. The cast iron, granular/polymer, crushable foam, MohrCoulomb, DruckerPrager/Cap, and jointed material models provide nonassociated flow with respect to volumetric straining and equivalent pressure stress. Since the flow rule and the hardening evolution rules are written in rate form, they must be integrated. The only rate equations are the evolutionary rule for the hardening, the flow rule, and the strain rate decomposition. The simplest operator that provides unconditional stability for integration of rate equations is the backward Euler method: applying this method to the flow gives AcpV= Y AA, agi A og (3.16) and applying it to the hardening evolution equations gives AHa = AA ,hi,a (3.17) where hi,, is the hardening law for H, . The strain rate decomposition is integrated over a time increment as: A = Acel + A (3.18) where As is defined by the central difference operator: S= sym Ax (3.19) a \x, + AAx The total values of each strain measure as the sum of the value of that strain at the start of the increment, rotated to account for rigid body motion during the increment, and the strain increment, are integrated. This integration allows the strain rate decomposition to be integrated into: 8 = C el+ pl (3.20) From a computational viewpoint the problem is now algebraic: the integrated equations of the constitutive model for the state at the end of the increment, must be solved. The set of equations that define the algebraic problem are the strain decomposition, the elasticity, the integrated flow rule, the integrated hardening laws, and for rate independent models, the yield constraints: fi = 0 (3.21) For some plasticity models the algebraic problem can be solved in closed form. For other models it is possible to reduce the problem to a one variable or a two variable problem that can then be solved to give the entire solution. For example, the Mises yield surfacewhich is generally used for isotropic metals, together with linear, isotropic elasticityis a case for which the integrated problem can be solved exactly or in one variable. 54 For other rateindependent models with a single yield system the algebraic problem is considered to be a problem in the components of AsE'. Once these have been found the elasticitytogether with the integrated strain rate decompositiondefine the stress. The flow rule then defines AA and the hardening laws define the increments in the hardening variables. CHAPTER 4 CONSTITUTIVE MATERIAL LAWS For bulk materials a constitutive highly nonlinear material model needs to be used that describes the solidlike behavior of the material during small deformation rates as well as the fluidlike behavior during flow conditions. Most of the models describing bulk solids were developed in civil engineering to describe sand, clay and rock behavior. Due to difference in stress magnitudes, loading and flow, only some of these models can be applied to bulk solid applications. They are based on a consistent continuum mechanics approach. The simplest model to use would be an elastic model based on Hooke's law. This model was used by Ooi and Rotter (1990). Other elastic concepts such as nonlinear elasticity (Bishara, Ayoub and Mahdy, 1983) and hypoelasticity (Weidner, 1990) have also been used. However these models cannot predict phenomena such as cohesion, dilatancy, plasticity and other effects that are essential for the model. These effects need to be described by elasticplastic models that can describe irreversible deformations. One of the first plastic models used is the MohrCoulomb model (1773). S1I= C2= C3 MohrCoulomb Tresca CY C3 Tr Figure 4.1: Tresca and MohrCoulomb yield surfaces It is derived from the Tresca criteria for metal plasticity and states that a material can sustain a maximum shear stress for a normal load in a plane. If the shear load is more than that maximum stress then the material starts to flow. Figure 4.1 shows the yield surfaces for both MohrCoulomb and Tresca in the principle stress space. Another yield surface that is numerically easier to handle is the DruckerPrager (1952) surface because it does not have the edges of the MohrCoulomb. Because of this, researchers have a preference to use this model. It is derived from the vonMises criteria and is shown in Figure 4.2 C 1= C2= 3 C72 DruckerPrager \ Von Mises C03 Figure 4.2: Von Mises and DruckerPrager yield surfaces One disadvantage to the above models is that under hydrostatic loading plastic deformation can not be predicted. That is because during isotropic loading (on the line CI= 02= 03), the stress lies within the yield surface. However, it is known from experiments that plastic deformation does occur during isotropic loading. Therefore a cap is introduced to create a closed yield surface that limits elastic deformation in hydrostatic loading. So, as can be seen from figure 4.3, the model consists of a yield cone and a yield cap. C1 02 3 02 C3 1 Figure 4.3: Closed yield surface In the above elasticplastic models, both the elastic behavior inside the yield surface and the plastic flow rule have to be determined. The flow rule can either be associated, where the shape of the yield surface determines that flow direction or non associated, where an additional surface, the flow potential, is defined. The general form of the flow rule is usually assumed to be potential where the strain increments are related to the stress increments by the following relationship: S= d G (a) (4.1) de" =dA ' where: dEs = plastic strain increase dA= proportionality factor G = plastic potential o, = Cauchy stress state The flow rule is associated if the potential G is equal to the yield function F. If they are not equal then it is nonassociated. Capped DruckerPrager Model (1952) The capped DruckerPrager model is used to model cohesive materials that exhibit pressuredependent yield. It is based on the addition of a cap yield surface to the Drucker Prager plasticity model, which provides an inelastic hardening mechanism to account for plastic compaction and helps to control volume dilatancy when the material yields in shear. It can be used in conjunction with an elastic material model and allows the material to harden or soften isotropically. It is shown in Figure 4.4. Plastic Transition surface, Ft Shear Failure, Fs LV Hardening e \ \ Hardening Elastic \ 11 Figure 4.4: The linear DruckerPrager cap model 3 3 +23 The DruckerPrager failure surface is written as: F = q p tan / d = 0 (4.2) where P and d represent the angle of friction of the material and its cohesion, respectively. p trace (u ) is the equivalent pressure stress 3 q J s is the misses equivalent stress r 9 S3 is the third stress invariant is the deviatoric stress S = c + PI is the deviatoric stress The cap yield surface has an elliptical shape with constant eccentricity. The cap surface hardens or softens as a function of the volumetric inelastic strain. It is defined by the following equation: 2 F= [P P]+2 + Rp R (d + p, tanp) = 0 (4.3) P FI+a a/ 8 where, R is a material parameter that controls the shape of the cap and a is a small parameter that controls the transition yield surface so that the model provides a smooth intersection between the cap and failure surfaces: 2 F [PPa]2+LpK (d+P ,tani) a(d+ptan)=0 0 (4.4) As can be seen from the equations above the model itself needs 7 parameters or measured material properties for calibration. In addition to describe the elastic behavior two parameters are needed (E Young's modulus, and vpoisson's ratio). The coefficient of wall friction [t also needs to be determined. 3D General Elastic/Viscoplastic Model In order to be able to predict the mechanical behavior of cohesive powders in silos or hoppers, a physically adequate constitutive model has been developed (Cristescu, 1987). In general, the constitutive model captures all major features of material mechanical response and accurately describes the evolution of deformation and volume change. Dilatancy and the related effects, such as microcrack formation damage, and creep failure can be predicted. The threedimensional general elastic/viscoplastic constitutive equation with non associated flow rule (Cristescu, 1987) is i 1 1 1 W(t) 8F E= +  l+k, 1 (4.5) 2G 3K 2G) H(, Tr) 8(5) = rate of deformation tensor H(cr, r) = yield function a = Cauchy stress tensor W' (t) = irreversible stress work per unit volume, a = mean stress F = viscoplastic potential, K, G = shear and bulk moduli k, = viscosity coefficient for transient creep W1 (t) 1 is the identity matrix. The function 1 W') is chosen to represent mechanical H(, z) behavior of the material due to transient creep and symbol ( ) is known as Macaulay bracket: (A)=(A+ JAI). The last term in equation 4.5 describes the mechanical behavior of the elastic/viscoplastic material that exhibits viscous properties in the plastic region only. In the elastic/viscoplastic formulations stress and strain are timedependent variables, and time is considered as an independent variable. In this problem the phenomena of creep and plasticity cannot be treated separately as only the superposition effect. It is worthy to note that the model initially included a term for steadystate creep given by the equation: . C 9S s = ks (4.6) where S(oa) is the viscoplastic potential for steadystate creep and ks is a viscosity coefficient. This term however describes behavior over long periods of time such as creep in tunnels. However, transient creep can describe the behavior of powders in hoppers at relatively smaller stresses and hence the above term is ignored. The material is assumed to be isotropic and homogeneous. Hence the constitutive equation depends on the stress invariants; primarily the mean stress: 1 ( =1 +02 + +03 (47) and the equivalent stress (T or the octahedral shear stress Z : 2 2 2 2 T=U = _1 +_2 +3 1 2 1 3 2 3 (4.8) or in terms of an arbitrary (1,2,3) coordinate system: = (11 22)2 +(11 33) +(02 033 2 6 22+ 23 +6132 (4.9) The strain rate consists of elastic and irreversible parts: = E + (4.10) The elastic part of the strain rate tensor is S= +  1 (4.11) 2G 3K 2G) (4.11) and the irreversible part is S= k 1 ) O(4. 12) Equation 4.11 describes an instantaneous response of the material to loading. The irreversible part, (equation 4.12), depends on the time history as well as on the path history of the stress. One of the time effects that the model describes is creep. Deformation due to transient creep stops when the stabilization boundary is reached. The equation of this boundary (the locus of the stress states at the end of transient creep) is H(o, ) =W'(t) (4.13) The irreversible stress work per unit volume W'(t) describes the irreversible isotropic hardening: W (t) = r:dE' (4.14) One of the concepts described by the model is whether the bulk solid is compressing or dilating. Figure 4.5 shows these domains. Triaxial experiments on cohesive materials show that powders show compressibility at lower stresses followed by dilatancy. Between the two domains we have the compressibility/dilatancy boundary. In the case of a triaxial test the boundary is determined by determining the where the slopes of the plots of the stress vs. volumetric strain are vertical. FAILURE DILATANT I CO 4P E 1S3B E COMPRESSIBLE Figure 4.5: Domains of compressibility and dilatancy G The viscoplastic flow occurs when: K> W'(t) > 0 (4.15) This model will enable us to capture all major features of powder mechanical response and accurately describe the evolution of deformation and volume change. It can predict the dilatancy and the related effects, such as creep failure. It is calibrated by numerous hydrostatic and deviatoric triaxial testing. Numerical Integration of the Elastic/Viscoplastic Equation The main difficulties in implementation of the elastic/viscoplastic constitutive model in the finite element code is the numerical integration of the irreversible (viscoplastic) part of the strain tensor. Suppose that at the time moment t all values of stress and strains in 4.5 are known. Consider the moment of time t + At. First, we represent the increment of the irreversible strain rate in the form of the truncated Taylor series with respect to the stress and irreversible strain increments: Ai' = (t)a + {(t)AE', (4.16) where (= W(t) )1F (4.17) SH(c, T) /O and W' (t + At) = W (t) + (t) ( (t + At) 1 (t)). (4.18) The total strain increment consists from elastic and viscoplastic parts. The elastic part is E =EE (4.19) and from the Hooke's law (==C E. ) (4.20) where C is constant stiffness matrix dependent on the elastic constants. The increment of irreversible strain is As = I (t + At) 7 (t) (4.21) and A' =^ ((1a) e(t)+c (t+)) (4.22) If a = 0, the explicit scheme (Euler forward scheme) for the integration of viscoplastic strains results. On the other hand, if a = 1, the fully implicit scheme (Euler backward scheme) for the integration is obtained. The case a = results in the socalled "implicit 2 trapezoidal" scheme. Substituting (4.22) in (4.21): = ij (t)+a (t)Aa + a (t)AEI (4.23) At 0G de or E (t)At + aAt (t) (a (t + At) a (t)) Ae' (4.24) 1 aAt (t)e t Returning to formula (4.20) the expression for stress tensor at the moment t + At gives: ,(t+AO) I {C(loN(t)((t+AO)c't))A 0 1 2V(t)+a (t) a1(4.25) The important issue to be addressed is stability of the numerical integration scheme. Note that the Euler forward method does not give an accurate numerical 1 solution. For a < , the integration process can proceed only for values of At less than some critical value (this critical value should be determined in some way). The accuracy of the numerical integration scheme with respect to stability and convergence has been estimated by comparing the theoretically predicted values for irreversible stress work and irreversible strain rate (obtained analytically by direct integration of the constitutive equation) with those obtained numerically during creep. There is good correspondence of results as for the irreversible stress work as for the irreversible strain rate. The elastic component is calculated from Hooke's law: o, = A25 + 2/E1 (4.26a) or a = C (4.26b) where 1 l, i = j, = 1= (4.27) i and pu are Lame constants. In another form 66 oa A +2p A A 0 0 0 ~, Y l A A+2pl A 0 0 0 yy Cz A A A+ 2p 0 0 0 (4.28) CXY 0 0 0 2p 0 0 exy oxz 0 0 0 0 2p 0 E5 CY 0 0 0 0 0 2p _L s Inverse matrix C' may be calculated using using the C matrix above. 8 = C1 ( (4.29) Relations of Lame's constants the bulk and the shear moduli are 2 K = i + u (4.30) 3 E= (3 ) (4.31) G = u (4.32) CHAPTER 5 EXPERIMENTAL RESULTS AND DISCUSSION Parameter DeterminationShear Tests Both direct shear testers and indirect shear testers are used to determine the parameters that are needed for the model per the procedures described in Chapter 2 (Saada, 2005). Silica powder with particle size of around 50 microns was used. To determine the internal angle of friction, /. and the cohesion, d the linear failure surface was determined using several Schulze tests. The Schulze test (1994) was chosen over the Jenike test (1961) because of the simplicity of the sample preparation and the test procedure. The consolidation weights ranged from 2kg to 25kg. The specimen was sheared at values between 20% and 80% of YeldLout Figure 5.1: Output data from Schulze test on Silica at 8Kg the steadystate value. Figure 5.1 shows the results of a consolidation weight of 8Kg. The slope of the yield locus (straight red line) is the internal angle of friction and the Y intercept is the cohesion. Plots of cohesion and internal angle of friction are plotted vs the principal stress o in Figure 5.2. Table 5.1 summarizes the results of 16 tests (Saada, 2005). The tests were also repeated to check for consistency. Table 5.1: Summary of Schulze test results i g ka3 (internal Weight (g) oa (kPa) intn) d (cohesion) friction) 2000 2000 3000 3000 4000 4000 6000 6000 8000 8000 12000 12000 16000 16000 22000 25000 1.57 1.55 2.29 2.36 3.05 3.13 4.55 4.64 6.02 6.50 9.16 9.08 11.97 11.98 15.99 18.24 28.71 27.87 29.00 28.80 31.29 30.37 34.86 33.18 31.47 30.07 32.33 33.83 34.09 33.65 32.18 32.05 0.359 0.364 0.43 0.5 0.47 0.566 0.46 0.639 0.75 0.849 1.04 0.935 1.11 1.16 1.38 1.56 The cohesion ranged from about 0.36KPa to about 1.56Kpa at a o0 value of 12Kpa. The data can be fit with a second order polynomial as seen in Figure 5.2 d (Cohesion) p(internal friction) 30 ,20 0 0 SOo 5 a,(Kpa) 10 ai (Kpa) Figure 5.2: Schulze test results for to determine the linear DruckerPrager surface 1.4 1.2 1 f0.8 _n * x x a 69 The internal angle of friction is relatively constant around the 32 deg value except foro, values that are less than 3Kpa where the angle drops slightly (Figure 5.2). To determine the wall friction coefficient, Jenike tests are preformed between 304 stainless steel slab, 120mm x 120mm (with 2B finish) and silica (Saada, 2005). The test was conducted using the procedure described above with weights decreasing from 16kg down to 2kg. To check for repeatability the test was repeated twice. Results are shown in 30  29  M 28  27  0 26   Test 1 25 ....... .......Test 2 0 24 u 23 D 22 0. .. .... 21 20 0.0 5.0 10.0 15.0 20.0 25.0 Normal Stress (kPa) Figure 5.3: Jenike wall friction results to determine boundary conditions Figure 5.3. It can be seen that the values of wall friction angles are slightly higher under the 5 kpa, but decrease in value and reach a constant value of about 22 deg as the normal stress is increased. This provides a value of u =0.4 to be used in the coulomb criteria. In the stress area of interest the test are repeatable and the maximum error is less than 5%. Detailed studies of wall friction are outside the scope of this research but are currently undertaken by other research groups. Future work in the area of simulation involves the study of more complex frictional models. Triaxial tests are required to determine the parameters for the cap in the model and the elastic parameters. A hydrostatic test where the sample is pressurized in all directions is preformed and the pressure and volume changes are recorded. Several deviatoric tests covering the range of confining pressures of interest are also needed. In these tests a fixed confining pressure is applied while the sample is given axial deformation. The stress, axial and volumetric strains are measured. Figure 5.4 shows results for five hydrostatic tests up to 69.0 kPa (Saada, 2005). The confining pressure is increased at a constant rate up to five values: 6.9KPa, 34.5KPa, 41.4KPa, 55.2KPa, and 69.2KPa. It is important to 70 60 0. 40 4 000 0 o)K 6.9 Kpa u030 Xj.rb A 34.5 Kpa )b so@.a0 41.4 Kpa 20 #0z x 55.2 Kpa 10Y o 69.0 Kpa 0 0.01 0.02 0.03 0.04 0.05 0.06 Ev Figure 5.4: Hydrostatic triaxial testing on Silica Powder (5 confining pressures) note that the specimen had an initial axial stress of 3.7KPa prior to the initiation of hydrostatic testing due to the filling procedure described in Appendix A. This is reflected in Figure 5.4 where the plots of the stressstrain curves show a clear jump at the beginning of the hydrostatic testing. The hydrostatic tests are followed by deviatoric tests where the confining pressure is kept constant and the axial deformation is increasing. 71 Figures 5.5 and 5.6 show the axial strain and volumetric strain respectively of five deviatoric tests. 250  oCa3(kpa) 200 AAAA 150  10 A 00000 A ^ 0000000 100 0, 000000000000 8 0o0o00 x34.5Kpa 50 A 4 o 13.8Kpa & 46.9Kpa Or IU o 4.8Kpa 0.005 0.045 0.095 0.145 0.195 0.245 8i Figure 5.5: Axial deformation of 4 deviatoric triaxial tests on Silica Powder and a rate of 0.1 mm/min 250  o1G3(kpa) 200 A A A A 150 Xx o x 34.5Kpa 100 D 0 Ao^*13.8Kpa 06, 53 4.8Kpa 0 50 X o o oY00 A46.9Kpa 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 ev Figure 5.6: Volumetric deformation of 4 deviatoric triaxial tests on Silica Powder and a rate 0.1mm/min As expected the graphs show an increase in measured axial stress with increasing confining pressure. Each test shows an increase in strength until a plateau is reached. The volumetric strain graphs show a compaction region followed by a dilation region. The data has been corrected for the confining effects of the latex rubber membrane using Hooke's law and the elastic properties of the rubber. To check for repeatability several tests were repeated more than once. Figure 5.7 shows the axial strain results of two deviatoric tests at 13.8KPa and 34.5Kpa. 250  oao3(kpa) 200 150  Figure 5.7: Axial deformation of 4 deviatoric tests on Silica Powder at 13.8KPa and 34.5Kpa50 34.5Kpa has is less than 50. 13.8Kpa 0.005 0.045 0.095 0.145 0.195 0.245 Figure 5.7: Axial deformation of 4 deviatoric tests on Silica Powder at 13.8KPa and 34.5Kpa It can be seen that the error in the shear stress is maximum at the higher deformations and has is less than 5%. Because of the low pressures that are being investigated, the plots show results that require the 'modification' of the testing procedure. This is because of the apparent oscillations in the graphs and the not very distinct value of failure for the specimen. Another issue to resolve is the sharp jumps and drops in stresses under 10kpa (shown in the circle), which are due to the loading piston friction while it is lining up with the center of the specimen. These oscillations are described and explained in Appendix A under piston friction error. In spite of these challenges these results do produce values that are satisfactory for use in the model. Diamondback hopper measurementsTekscan Pads The experimental procedure for the measurements of wall stresses in the hopper was described in chapter 2. Measurements were taken when the powder was static and at belt speeds of 40, 60, 70, and 100 (Johanson, 2003). Figure 5.8 shows the geometry of the mechanism at the outlet of the hopper . Powder in chute Pow der on belt B l S; Belt cr Figure 5.8: Discharge mechanism and belt Based on the dimensions of all the parts and the outlet, the velocity of the belt, table 5.2 was created to relate velocity and flow rate of the powder in the chute. Table 5.2: Belt velocity vs. flow rate Belt velocity Flow rate (cm3 /sec) 40 35.6 60 52.7 70 60.9 100 80.6 Figures 5.9 shows the location of the pad in the hopper and an output from a test. The area that is covered by the rectangles represents the line in the hopper where the transition between the vertical and the converging parts of the hopper. The left side of the pad starts midway through the flat section and loops around the curved section. Figure 5.9: Pad location and output example Rectangles are used to obtain wall stress as opposed to lines so as to obtain an average and not localized stresses. The data were averaged along eight paths shown by the white lines in figure 5.10 (Saada, 2005). They are 2.8cm apart and span from the middle of the flat section of the hopper to halfway around the curved section (90 deg). Figure 5.10 shows the contact pressure measured while the powder was static at a depth of 2.8cm from the hopper transition. These stresses generally decrease around the perimeter and are smaller at greater material depths. The results show a peak at the flat section of the hopper and oscillations going around the curved section of the hopper. I Wall Stress (@2.8 cm) Static Initial 2 cm ' ^ LgA M~ LI~t 40 Angle (deg) Figure 5.10: Static Wall stress data at a distance 2.8cm below hopper transition Figure 5.11 shows wall stress measurements were along the same path of depth 2.8cm conducted at various speeds 60, 70 and 100 in addition to the static state (Saada, 2005). The plots show agreement in the pattern and magnitude of wall stress at the various belt speeds. 14 12 10 \ Initial M (M .speed 60 8 uspeed 70 6 speed 100 U,4 2 0 20 40 60 80 Angle (deg) Figure 5.11: Wall stress data at a distance 2.8cm below hopper transition at various speeds and static conditions Figure 5.12 shows the same data at the same speeds at a depth of 5.6cm. Appendix A shows the results at the rest of the six paths shown in Figure 4.10. 10 9 Initial 8 *speed 60 S7 \ speed 70 speed 100 (43 1 0 0 20 40 60 80 Angle (deg) Figure 5.12: Wall stress data at a distance 5.6cm below hopper transition at various speeds and static conditions Figure 5.13 shows the results of a Fourier series analysis conducted at the results from two different regions (Johanson, 2005). Region numbered #1 is in the flat section of the converging hopper. In this region the fluctuations in the stress have a wide range of frequencies. Region #7 is in the round section. In this region the fluctuations have a wide range of frequencies. It is hypothesized that the converging diverging nature of the Diamondback hopper causes a shock zone formation across the hopper that propagates at points along the edge of the flat plate section. Here the radius of curvature changes causing a stress discontinuity. This should be a region of steep stress gradients. This is validated by experiment it will be compared to FEM calculation. 77 0.(0 Radius of curvature 0.6 0.5 changes causing a 0.5 stress discontinuity. R o0.(4 This is a region of O.o 4 steep stress B 0.0 5. Ogradients. . 10 100 11 10 100 i. Period of Signal (sec) Period of Signal (sec) 0 500 1000 1500 .. . .Time (se 0 500 1000 1500 2001 Time (sec) Time (sec) Figure 5.13: Fourier series analysis on the wall pressure measurements. Conclusions and Discussion of the Test results and Testing Procedure Both simple direct shear tests and more complex triaxial tests are needed to determine the parameters for the models that are used in this research. Sixteen Schulze tests were conducted to determine the internal angle of friction, 3, and the cohesion, c (or d). The tests were conducted at various consolidation stresses and were repeated at least twice to check for repeatability. The internal angle of friction was determined to be constant at a value of 320 and the cohesion was increasing quadraticaly. Jenike tests were conducted to determine the wall friction angle. That angle was determined to be 220. Both hydrostastic tests and deviatoric triaxial tests were performed for the calibration of the viscoplastic model. The deviatoric tests were conducted at a confining pressure of up to 47Kpa. They were repeatable and showed the expected increase in strength with increasing confining pressure and the showed the powder going from compressibility to dilatancy. Performing the tests at this lower confining pressure (less than 50Kpa) created issues that could introduce some error to the results. This is partially seen by the fact that the material does not show an apparent failure point. The two main issues are: 1) the influence of the rubber membrane on the deformation of the specimen, and 2) the influence on the confining pressure of the water height. This 2nd issue is a problem due to the fact that the stress (due to the water pressure) at the top of the sample is close to zero while at the bottom it is around 2KPa which could present some error especially at the test with 4.7Kpa confining pressure. For the 1st issue, Appendix B provides a procedure to correct for membrane pressure based on hoop's law and the rubber properties. To minimize the influence of the membrane it is suggested that a much thinner membrane is used. This, however, may cause some complications in the 'current sample preparation procedure', since extra care has to be taken so as not to cause any damage to the membrane. It is also suggested that plastic wrap be used instead of rubber membrane, since it is thinner and because of plasticity, it would be easier to correct for it. The ideal procedure would be to conduct the test without a membrane but since this is impossible with the current setup it is suggested to use another tester such the Cubical Biaxial test. CHAPTER 6 FEM RESULTS AND DISCUSSIONS Capped DruckerPrager Model The capped DruckerPrager model is used to model cohesive materials that exhibit pressuredependent yield. It is based on the addition of a cap yield surface to the Drucker Prager plasticity model, which provides an inelastic hardening mechanism to account for plastic compaction and helps to control volume dilatancy when the material yields in shear. It can be used in conjunction with an elastic material model and allows the material to harden or soften isotropically. It is shown in Figure 6.1. q o3 Plastic Transition surface, Ft Shear Failure, Fs \ 2 \Cap, F 1 \\ Hardening 1 \\ \ Hardening 'o Elastic \ P (, + 203) Figure 6.1: The Linear DruckerPrager Cap model The DruckerPrager failure surface is written as: F = q p tan 8 d = 0 (6.1) where 3 and d represent the angle of friction of the material and its cohesion, respectively. p L Ltrace (a ) is the equivalent pressure stress 3 q is the misses equivalent stress 2 S = a + PI is the deviatoric stress The cap yield surface has an elliptical shape with constant eccentricity. The cap surface hardens or softens as a function of the volumetric inelastic strain. It is defined by the following equation: 12 F,= 7[P Pa 2 + Rp R (d + p, tan/ )= 0 (6.2) S / os 8 where, R is a material parameter that controls the shape of the cap and a is a small parameter that controls the transition yield surface so that the model provides a smooth intersection between the cap and failure surfaces 2 f= [pp]2+[p1 (d+p.tanf) a(d+p tan/) =0 cosa p1 (6.3) As can be seen from the equations above the model itself needs 7 parameters for calibration. In addition to describe the elastic behavior two parameters are needed (E Young's modulus, and vpoisson's ratio). The coefficient of wall friction [t also needs to be determined. Results from the shear testers listed in Chapter 5 are used to determine the required parameters and are listed in Table 6.1. Table 6.1: Capped DruckerPrager parameters determined from shear tests E u p c (or d) R a K Evo Pa ev 4.4 0 Run 1 4.6x107 0.24 33.7 1.16 0.33 0.04 1 0.0024 10 0.0086 20 0.021 30 0.0321 0.5 0 Run 2 4.6x107 0.24 28.7 0.359 0.33 0.04 1 0.0001 10 0.0086 20 0.021 30 0.0321 Drucker Prager Model Verification The DruckerPrager model was used to perform FEM simulations on the results of a triaxial test. The results are used to compare the way perfect plasticity approximates powder behavior. The specimen was simulated as a deformable body (height=15cm, diameter=7cm) and discritized into 1063 tetrahedral elements as seen in Figure 6.2 a). A confining pressure of 10 OKPa was applied around the specimen and the specimen was deformed axially by 3cm (up to 20% strain) as seen in Figure 6.2 b). a) Undeformed snecimen b) Deformed snecimen Figure 6.2: FEM simulation of a triaxial test using the DruckerPrager Model The parameters for the model where determined from the shear tests and are shown in Table 6.1. The axial stress was calculated from the average of the stress at the top of the specimen. The axial stressstrain curves from the simulation were compared with the results that were obtained experimentally as seen in figure 6.3. 40 Silica @1lmnvmin 35  30 25  20 ** ***** * 15 Experimental 5 0 0 0.02 0.04 0.06 0.08 0.1 Strain Figure 6.3: Axial stressstrain curves for experiment and FEM using DruckerPrager the model Drucker Prager Model Hopper FEM Results The results from the shear tests above were used to run an ABAQUS FEM simulation on the diamondback hopper in the static case (during storage). As described in section 3.2, the powder is simulated as a deformable body with quadratic tetrahedron (C3D10OM) elements. The powder was divided into 11 different sets and the gravitational force was applied to each set to simulate filing. The steel wall is simulated as a discrete rigid body that required discritization. The boundary conditions where the tangential Coulomb friction of [t=0.4 on the walls and a vertical displacement of zero at the outlet. The loading is a gravitational constant of 9.81m/s2. Figure 6.4 shows a side view of the calculated contact stresses at the end of each filling step (Saada, 2005). They show the increase in stress at the bottom of the hopper and moving up as the force is applied. Figure 6.4: FEM calculated contact stresses at various filling steps It is important to note that the stress concentrations in the vertical section of the hopper (above the transition to a converging section), are not a powder phenomena but are rather due to FEM error in the contact analysis procedure. This is due to the fact that some of the elements in the powder penetrate the wall and cause stress concentrations. Also the nonsymmetry in the simulation could be due the fact that the meshing is not 100% symmetrical (automatic meshing techniques were used). Another explanation is probably due to the fact that the simulation of filling is not completely static and there is movement of powder nodes on wall nodes causing nonsymmetry. To simulate discharge the surface at the bottom of the powder was moved by a distance of 5cm. Figure 6.5 shows the contact stresses at 4 different intervals (Saada, 2005). Note Figure 6.5: FEM calculated contact stresses at various discharge steps the stress change close to the outlet and in the converging section. The comments above explain the nonsymmetry of the simulations. The area where the stress measurements are being compared and analyzed is shown in Figure 6.6. The FEM section shown shows the nodes where contact stresses are recorded. They are then compared with the Tekscan measurements along the corresponding paths shown as white lines in the plot. TekScan FEM Fi ...'^ t.. ... .* *. '.^ . Ir I. Figure 6.6: FEM vs. TekScan area of interest It is important to note that the reasons the white lines are curved is due to the fact that the pad is curved inside the hopper and is presented a flat section in the figure. Also due to the complexity of the contact analysis simulation, a balance between mesh density, element types and computer simulation time an appropriate mesh had to be chosen. The mesh chosen had 14 nodes across each of the paths that were analysed. These are seen as the black dots in the FEM contour plot. Figure 6.7 shows the comparisons of the static wall stresses at seven locations (from 2.8cm to 19.6cm below the hopper transition). The plots here show a similar trend. They capture the general magnitude of the wall stresses that are higher at the middle of the flat section but are decreasing as we go round in the curved section. The simulations show some oscillations, but fail to calculate the same oscillations that are recorded by the tekscan pads. The magnitudes of the stresses are better approximated at the flat section (between 0 and 40 deg) but are slightly higher at the range of 40 deg to 90 deg. It is possible that the oscillations in the pad readings are due to that fact that the powder is not completely static and there is some movement while the powder is settling and the measurement are being made. As explained earlier the simulation of flow was done by applying a displacement of 5 cm to the bottom surface of the powder. Another method attempted was the removal of the boundary conditions at the surface this however produced large deformations of the elements and caused the simulation to abort. The data analysis was conducted at various locations in the flow simulation. Figure 6.8 shows the comparisons at two seven paths SMeasured Static Wall Stress 0 FEM Static Wall Stress a 0 o I B w 45 Angle (deg) 90 M Measured Static Wall Stress 0 FEM Static Wall Stress 45 Angle (deg) 90 45 Angle (deg) 0 o o0 @14cm 45 Angle(deg) 90 0 45 Angle (deg) 90 45 Angle (deg) 90 000 0o @19.6cm A0 90 45 Angle (deg) 9 Figure 6.7: FEM calculated static wall stresses vs. measured wall stresses at seven locations Angle (deg) (@8.4 cm) OO 000 0 0 0 NIIII0 0 SMeasured Flow Wall Stress  X FEM Flow Wall Stress at 0.25 cm ........... ::G::::::: FEM Flow Wall Stress at 5 cm "R" 6 =5( S4) 3 45 Angle (deg) 90 Measured Flow Wall Stress  X FEM Flow Wall Stress at 0.25 cm O ii::::*: FEM Flow Wall Stress at 5 cm 45 90 Angle (deg) 45 Angle (deg) 90 45 Angle (deg) 45 Angle (deg) 90 12 45 Angle (deg) 90 10 @19.6 cm 6 4O 2 0 45 Angle (deg) 90 Figure 6.8: FEM calculated contact stresses vs. measured wall stresses at seven locations during flow and at two simulation displacements (0.25cm and 5cm). All the results show that the DruckerPrager model approximates the wall stresses during flow better than during storage. The stresses follow the same pattern of decreasing as the angle is increased. The simulation captures a good estimate of magnitude and some of the oscillatory pattern. We see that the magnitude of the oscillations increases as the displacement is increased. This is probably due to the increase in contact wall friction criteria between the nodes at the wall and the nodes at the surface of the powder. Several other parameters that affect the model were also investigated. The effect of flow rate is shown in Figure 6.9. The 5 cm displacement was applied at times of 1 second and 0.1 seconds. The figure shows the stresses at 5.6 cm and 8.4 cm transitions. MEASURED Flow Wall Stress I MEASURED Flow Wall Stress 9 X FEM at Discharge time=0.1 sec X FEM at Discharge time=0.1 sec "iiiiiiiiiiiii EM at Discharge time=1 sec 7 0 FEM at Discharge time=1 sec 8 7 6 0 20 40 60 80 0 I Anle (deg) 0 20 An4ge (deg) 60 80 a) At 5.6cm from hopper transition b) At 8.4cm from hopper transition Figure 6.9: FEM calculated contact stresses vs. measured wall stresses at two locations during flow at two different times 