UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 MODELING GR ANULAR HOPPER DISC HARGE AND SEGREGATION FOR WET COHESIVE PARTICLES By ANSHU ANAND A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2009 1 PAGE 2 2009 Anshu Anand 2 PAGE 3 To m y Mom and Dad 3 PAGE 4 ACKNOWL EDGMENTS Working for this PhD thesis has been some of th e best years of my life. I owe this pleasure to a lot many people. First and foremost, I would like to thank my advisor Prof. Jennifer Sinclair Curtis for her constant support, encouragement and guidance. Secondly, I would like to thank my research collaborators: Prof. Carl Wassgren at Purdue University, Dr. Bruno Hancock and Dr. William Ketterhagen at Pfizer, Groton for their invaluable suggestions, insight and support. Without the assistance from all these people this thesis would not have been possible. I would like to thank people in my research group at UF for valuab le discussions ranging from particle technology to sports, music and m ovies. During the course of my PhD work, it was always easy to find solutions to the problems f aced by me, whether it be an experimental or computational problem, with the help of these people. I would also like to thank a ll my roommates in these 4 years for all the good times and relieving the stress of work. Fi nally, I would like to thank my mom and dad to whom my thesis is dedicated. 4 PAGE 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........7 LIST OF FIGURES.........................................................................................................................8 ABSTRACT...................................................................................................................................10 CHAPTER 1 DISCHARGE RATE A ND SEGREGATION.......................................................................13 Introduction................................................................................................................... ..........13 Discharge Rate................................................................................................................. .......14 Segregation.............................................................................................................................19 Cohesion.................................................................................................................................20 Objectives...............................................................................................................................21 2 COMPUTATIONAL MODEL...............................................................................................24 Introduction................................................................................................................... ..........24 DEM Simulations...................................................................................................................24 Contact Detection............................................................................................................25 Contact Model.................................................................................................................25 Computational Domain........................................................................................................... 28 Liquid Bridge Model..............................................................................................................28 Simulation Pa rameters.......................................................................................................... ..32 3 EXPERIMENTAL................................................................................................................. .42 Introduction................................................................................................................... ..........42 Motivation...............................................................................................................................42 Materials and Experimental setup.........................................................................................42 Procedure................................................................................................................................43 4 RESULTS AND DISCUSSION.............................................................................................50 Introduction................................................................................................................... ..........50 Discharge Rate for Cohesionless Particles.............................................................................50 Fill Height and Friction...................................................................................................50 ParticleParticle and Pa rticleWall Friction.....................................................................52 Coefficient of Restitution................................................................................................52 Hopper Width..................................................................................................................5 3 5 PAGE 6 Hopper Angle ..................................................................................................................53 Outlet Width................................................................................................................... .54 Particle Diameter.............................................................................................................54 Particle Size Distribution.................................................................................................55 Discharge Rate for Cohesive Particles...................................................................................55 Effect of Outlet Width and Bond Number......................................................................56 Effect of Liquid Content..................................................................................................58 Effect of Hopper Angle...................................................................................................60 Segregation of Cohesive Particles..........................................................................................60 Effect of Bond Number on Segregation..........................................................................61 Effect of Liquid Content on Segregation.........................................................................63 Experimental Validation........................................................................................................ .64 Discharge Rate................................................................................................................. 64 Angle of Repose..............................................................................................................65 Clump Size..................................................................................................................... .65 5 CONCLUSION AND RECOMMENDATIONS.................................................................100 Conclusions...........................................................................................................................100 Recommendations................................................................................................................ .103 APPENDIX A SAMPLE IMPUT FILE.....................................................................................105 APPENDIX B HOPPER DISCHARGE CODE.........................................................................107 LIST OF REFERENCES.............................................................................................................242 BIOGRAPHICAL SKETCH.......................................................................................................249 6 PAGE 7 LIST OF TABLES Table page 21 List of baseline parameters for study of cohesionless disharge rate..................................38 22 List of hopper and particle propertie s used in experimental correlations..........................39 23 List of baseline parameters for study of cohesive disharge rate........................................40 24 List of baseline parameters fo r study of cohesive segregation..........................................41 31 Types of experimental cohesive system............................................................................49 41 List of various frictional cases........................................................................................... 98 42 Comparison of experimental discharge rate to that of simulations in the 90 hopper.......98 43 Comparison of experimental discharge rate to that of simulations in the 55 hopper.......99 44 Comparison of experimental results for angle of repos e to that of simulations................99 7 PAGE 8 LIST OF FI GURES Figure page 11 Outline of computational work..........................................................................................23 21 A schematic of normal and tangential contact force models.............................................35 22 Schematic of the quasi3D computational dom ain that models a thin slice of a large, rectangular crosssectioned hopper....................................................................................36 23 Representation of the liquid bridge fo rmed between particles of unequal size.................37 24 Illustration of the liquid br idge force as a function of th e separation distance between two particles.......................................................................................................................37 31 Schematic of the experimental hopper. A) Assembly B) Crosssection of the 55 hopper C) Crosssection of the 90 hopper D) Gate and orifice of the hopper..................45 32 Snapshots for calculating angle of repose. A) Experiments B) Simulation.......................48 33 Example of a snapshot from high speed camera for calculation of clump size.................49 41 Discharge profile for differ ent friction coefficients. A) = 0.05 and w = 0.05 and B) = 0.84 and w = 0.84. The value in parentheses is the mass discharge rate...................67 42 Instantaneous mass discharge rate pl otted as a function of instantaneous dimensionless fill height for different friction coefficients. A) = w = 0.05.B) = w = 0.1 C) = w = 0.15 D) = w = 0.2 E) = w = 0.5 F) = w = 0.84..................69 43 Mass discharged from the hopper plotted against time for a range of particleparticle and particlewall friction values. The value in parentheses is the mass discharge rate.....75 44 Mass discharged from the hopper plotte d against time for two values of the coefficient of restitution..................................................................................................... 76 45 Mass discharged from the hopper plotted ag ainst time for hoppers with two different widths, W ............................................................................................................................77 46 DEM mass discharge rate plotted against hopper angle (from the vertical). The correlations by Brown and Richards [Eq. 14] and Rose and Tanaka [Eq. 15] are also plotted................................................................................................................... ......78 47 Mass discharge rate per unit length raised to the 2/3 power plotted as a function of the outlet width, W0. The solid line is the trend predicted by the Beverloo correlation..........................................................................................................................79 8 PAGE 9 48 Mass discharge from the hopper plotted against time for two different particle diameters, d ........................................................................................................................80 49 DEM mass discharge rate plotted against the fines mass fraction for a bidisperse system. The correlation of Humby et al [Eq.17] is also shown......................................81 410 Dimensionless mass discharge rate as a function of Bond number for various values of dimensionless outlet width and L = 2%.........................................................................82 411 Mass discharge rate per unit length raised to the 2/3 power plotted as a function of the outlet width, W0, for various values of Bond number and L = 2 %.............................83 412 Beverloo constant, k as a function of Bond number for L = 2%.......................................84 413 Snapshot of discharge. A) Cohesionless, Bo = 0, Wo/d =10.6, L = 2% B) Cohesive, Bo = 2.5, Wo/d =10.6, L = 2 %..........................................................................................85 414 A comparison of DEM simulation results with predictions from the modified Beverloo correlation for cohesion......................................................................................87 415 Mass discharged from the hopper as a function for time for different Bond Number. The values in parentheses show the mass discharge rate per unit length. A) Bo = 0.3 B) Bo = 2............................................................................................................................88 416 Mass discharge rate plotted against hoppe r angle (from the vert ical) for different Bond numbers with Wo/d = 10.6 and L = 2 %...................................................................90 417 Discharge segregation for varying Bond numbers. L = 2%, f x = 5%, D = 4.3, and = 90 A) for Bo < 1 B) for Bo > 1................................................................................91 418 The normalized segregation index plotted as a function of Bond number........................93 419 Local fines fraction fields. A) Bo = 0 (c ohesionless), and B) Bo = 1. Red indicates a large fines fraction as compared to a well mixed state, and blue indicates a small fines fraction. Gray regions indicate well mixed regions.................................................94 420 Discharge segregation for varying liquid content (L), with Bo = 0.5, f x = 5%, D = 4.3, and = 90.................................................................................................................95 421 Illustration of Clump sizes from simulation. A) Bo=0.3 B) Bo=1 C) Bo=3.5...................96 422 Illustration of Clump sizes from expe riments. A) Silicone oil with Cast iron, Bo=0.3 B) Water with Cast iron, Bo= 1 C) Silicone Oil with Glass, Bo=1 D) Water with Glass, Bo=3.5.....................................................................................................................97 9 PAGE 10 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy MODELING HOPPER DISCHARGE US ING DISCRETE ELEMENT METHOD By Anshu Anand August, 2009 Chair: Jennifer Sinclair Curtis Major: Chemical Engineering Science The present work investigates hopper discharg e phenomenon from a quasi3D, rectangular hopper using the Discrete Element Method. Accura te prediction of the discharge rate from hoppers is important in many industr ial processes involving the hand ling of granular materials. For cohesionless granular media, the effects of particle properties (particle size and size distribution) and hopper geometry (hopper width, outlet width, angle and fill height) are studied and compared to previously published experimental correlations. The resu lts indicate that DEM simulations are fully capable of reproducing trends in the disc harge rate that are wellknown experimentally. For example, particle size and hopper width are shown to have a minimal influence on the discharge rate. In addition, for rectangular hoppers, the di scharge rate is shown to vary with the outlet width raised to the 3/2 power as given by the modified Beverloo correlation. The DEM simulations are also used to explore a wider range of parameters that have not been or are not easily explor ed experimentally. For exampl e, the effects of hopper friction, particle friction, and coefficient of restitution are invest igated, and particle friction is shown to have a significant influence on the hopper discharge behavior. The present work also investigates the para meters affecting the discharge rate of a wet cohesive system. The cohesion between the partic les is described by a pe ndular liquid bridge 10 PAGE 11 force m odel and the strength of the cohesive bond is characterized by a Bond number. The Beverloo correlation is applied to cohesive sy stems by modifying the Beverloo constant as a function of Bond number. The predictions obtai ned from this modified correlation fit the simulation data reasonably well. In addition, the effect of hopper angle in cohesive systems is shown to follow a trend similar to cohesionless systems, where the discharge rate is insensitive to changes in hopper angle except be low a critical angle (with resp ect to the vertical) where the discharge rate increases rapidly. This critical angle of flow de creases with increasing cohesion. Granular materials may readily segregate due to differences in particle properties such as size, shape, and density. Se gregation is common in industr ial processes involving granular materials and can occur even after a material ha s been uniformly blended. One of the specific objectives of this work includes investigating via simulation the effect of particle cohesion due to liquid bridging on particle segregation. Specifica lly, a bidisperse granular material flowing from a 3D hopper is simulated using the discrete element method (DEM) for cohesive particles and the extent of discharge segr egation is characterized over time The strength of the cohesive bond is characterized by the Bond number determin ed with respect to the smaller particle species. As the Bond number of the system increase s, the extent of discha rge segregation in the system decreases. A critical value of Bo = 1 is identified as the condition where the primary mechanism of segregation in the cohesionless hopper system, i.e. gravityinduced percolation, is essentially eliminated due to the liquid bridges between particles. Finally, experiments are performed in a hopper of identical dimensions as in the simulation to verify the model. Three specific characteristic of the flow: discharge rate, angle of repose of the material and size of clumps formed during th e discharge are investigated experimentally and compared to the simulations. The simulation results agree well with the e xperimental results. In 11 PAGE 12 12 general, the discharge rate decr eases, angle of repose increases and size of clumps increases as cohesion increases. PAGE 13 CHAP TER 1 DISCHARGE RATE AND SEGREGATION Introduction Hoppers are important devices that are widely used during the processing and handling of granular materials. One of the major industrial problems in using hoppers is obtaining reproducible and consistent partic le discharge flow rates. This problem is associated with the complex flow patterns of granular materials insi de the hopper. Accurate prediction of discharge rate is necessary for dependable design and working of hoppers. For example, in the pharmaceutical industry, exact dosing and metering of the quantity of material flowing out of a hopper is critical in a die filling process for tablet production. Granular materials may readily segregate due to differences in particle properties such as size, shape, and density. Se gregation is a natura lly occurring phenomenon common in industrial processes involving granular mate rials and can occur even after a material has been uniformly blended. In most instances, this segregation is problematic in that product quality is usually dependent on maintaining homogeneity of th e blend. For example, in the pharmaceutical industry, segregation during the production of so lid dosage forms may cause variations in the active drug ingredient concentration from tablet to tablet. A deeper understanding of hopper discharge phenomenon is required to keep the above problems under control. The focus of this thesis work is to understand this phenomenon using the tools of computer simulations a nd experiments. The first part of this chapter discusses the discharge rate from hoppers followed by a discussion of granular segreg ation. Then, cohesion is discussed briefly. Finally, the objective s of this work are outlined. 13 PAGE 14 Discharge Rate One of the first studies of granular flow s in a hopper was based on a continuum model by Janssen (1895). His analysis showed that above a certain height, the wei ght of the particle bed within a hopper is supported by frictional forces on the side walls. Thus, the pressure on the bottom of the hopper is independe nt of fillheight. Since the overburden pressure has a significant influence on the discharg e rate through the orifice, Jan ssens theory predicts that the hopper discharge rate becomes fillheight indepe ndent for bed heights above a critical value. Significant efforts have been made since then to predict and model the discharge of granular material from hoppers using experiments and, more recently, through computer simulations. Numerous correlations have been developed over the years to predict th e discharge rate of bulk particulate materials from hoppers. Taking into account the fillheight independence of the discharge rate, a dimensi onal analysis suggests that where Wis the mass discharge rate,1/25/2 0WgW is the characteristic bulk density of the system, g is the gravitational acceleration, and is the outlet width. However, many of the early experim entalists found that the discharge rate varies with th e outlet width raised to a power slightly greater than 2.5. The experimental work of Franklin and Johanson (1955) yielded a power of 2.93, Rausch (1948) found the power to be 2.8, and Brown and Richards (1960) observed a value of 3.1. It was Beverloo et al. (1961) who were able to resolv e the apparent discrepancies. 0WBeverloo et al. (1961) worked with a variety of sands having different mean particle sizes in the range of 93 to 300 microns, and seeds like linseed, spinach, waterc ress, rapeseed, kale, sweed and turnip which had a variety of shapes and sizes. Beverloo et al. (1961) plotted agains t and obtained a straight line with an intercept Z on the abscissa. The value of the intercept Z is correlated to the mean particle diameter d and a parameter k, where Z = kd. This 2/5W0W 14 PAGE 15 intercept was used to m odify the outlet width to a ccount for the effects of a layer of particles just near the edge of the orifice wh ere there is little or no flow. Th e resulting correlation, referred to now as the Beverloo correlation, is among the most c ited particle discharge rate correlations and was originally developed for a flatbottomed, cylindrical silo with a circular opening: 15 2(flatbottombOWCgWk 2) d (11) where k is the Beverloo constant, wh ich is a particle shapedependent constant. The constant C is found e mpirically to be in the range 0.55 < C < 0.65. The quantity b is the bulk density of the material after filling the hopper. Huntington and Rooney (1971) showed that different degrees of compaction during filling had littl e effect on the discharge rate. They suggested using the flowing bulk density at the hopper outlet, flow, rather than the filled bulk density. Typically, f low is determined by taking the ratio of the mass flow rate to the volumetric flow rate, the latter quantity measured by observing the change in the material head in the hopper: flowW Q (12) Using this flow density in the Beverloo formulation, the value of C lies in the narrower range of 0.575 < C < 0.595. It should be noted that a typical value of 0.58 is used in most applications of the Beverloo equation. According to Beverloo, the difference between the experimental and predicted disc harge rates is approximately 5 percent on average and 12.5 percent maximum. The Beverloo correlation has since been modi fied for various hopper geometries. Myers and Sellers (1971) extended the correlation to a flat bottome d (90 degree) rectangular hopper outlet of width W0 and length L and predicted the discharge rate as: 15 PAGE 16 13 21.03()()flatbottom flow OWg L k d W 2k d (13) Note that th e discharge rate varies as 3 2(OWkd ) in this case instead of 5 2(OWkd ) as for the case of hoppers with circular openings. It is dimensionally still consistent because of the ( Lkd ) term for the rectangular hoppers. To take into account the e ffect of the hopper angle, Brown and Richards (1965) analytically modified the co rrelations developed for flatbottomed hoppers as follow: 3 2 5 21cos sinflatbottomWW (14) where is the hopper angle from the vertical and f latbottomW is the discharge rate for a flat bottomed hopper with the same hopper width and outlet width. Rose and Tanaka (1959) proposed an alternate empirical method for taki ng into account hopper angle effects. They recommend the use of a multiplicative correction factor F for discharge from conical hoppers: 0.351 tantan1 (tantan) tantan1 F (15) where is the hopper angle fr om the vertical and is the angle from the horizontal of the stagnantflowing boundary, which is constant for a given particulate material. This angle is determined in situ by observation. It is important to note that this angl e is different from the angle of repose of the material. A value of = 45 is typically assumed in the absence of any other information. Both hopper angl e correlations predict an increa se in the discharge rate with decreasing hopper angle, with th e discharge rate increasing rapi dly as the hopper angle tends to zero. 16 PAGE 17 All of the previous correlations were devel oped f or flow of nearly monodisperse granular materials. Humby et al. (1998) analyzed the experimental work of Arteaga and Tuzun (1990) in which binary mixtures of acryl onitrile butadiene styrene (ABS), acrylic, radish and turnip particles were used (sizes ranging from 600 m to 3 mm) and modified the Beverloo correlation for a bidisperse system by including a parameter Z in place of the kd term in the Beverloo equation. Z depends on the volume, moment mean diameter. dVM that is a function of the number fraction xi of the large and small particles and the diameter di of the large and small particles ( i i d i x i i d i x VM d 3 4 ). Z is given as 0009.0 1939.2)(' VM d VM dZ (16) and Z decreases as the fraction of fines increas es. Taking into account the flowing density of the binary mixtures, the Beverloo correlation is modified as 13 21.03()()flatbottom flow OWg L k d W 2Z (17) where f low and Z are the modified parameters for the binary mixture. There is a slight increase in the discharge rate with increasing fines mass fraction which can be attributed to a closer packing of particles and an increased flowing density. Most researchers have noted only a slight dependence of the disc harge rate on the fill height, though Newton et al. (1945) reported that the disc harge rate is proportional to where H is the f ill height. Typical experimental results show that discharge rate is independent of fill height above a critical height. Myers and Sell ers (1971) suggested that the discharge rate is constant provided that fill height at the cen terline is greater than the orifice diameter W0. 0.04H 17 PAGE 18 There is also a general agreem ent that the discharge rate is independent of hopper width as long as it is not too sm all. Brown and Richar ds (1960) reported that the discharge rate is constant provided W > 2.5W0, where W is the hopper width and W0 is the outlet width. For smaller values of W larger discharge rates are observed and in the limit the whole mass accelerates indefinitely u nder grav ity (Nedderman et al. 1982). 0WWIn general, particle diam eter has been reported to have little effect on the discharge rate as long as the orifice is sufficiently large compared to the particle diameter. For orifices less than approximately six particle diameters, the flow is intermittent and irreproducible (Nedderman et al. 1982). The Beverloostyle correl ations for predictions of the discharge rate cannot be used in such cases. For spherical particles, the Beverl oo equation predicts that smaller particles will have a slightly larger discharge rate due to the smaller value of the correction factor kd Hence, the effective flowing discharge area for smaller part icles is larger and a larger discharge rate is expected; however, this effect te nds to be relatively minor. Researchers have reported that there is littl e to no effect of par ticlewall friction on the discharge rate. The effect of wall roughness was studied by Nyugen et al. (1979). They reported a slight decrease in the flow rate with increasing particlewall friction, but this friction rarely affects the discharge rate by more than 10%. There have been no comprehensive experiment al studies investiga ting the effect of particleparticle friction or co efficient of restitution on hopper di scharge rate. There is an implicit assumption in all of the experimental correlations that these interaction parameters play little or no role in determining the discharge rate. The previous correlations are not applicable for fine particulate materials in which interstitial air effects play a crucial role in de termining the discharge ra te. Several correlations 18 PAGE 19 exist in the lite rature to model such systems [(Crewdson et al. 1977), (Verghese and Nedderman, 1995) and (Spink and Nedderman, 1978)]. Interestingly, even in these correlations for fine particulates, the discharg e rate still varies with the outlet width raised to the power of 5/2. There are a few other correlations that exist in the literature [(Barletta et al. 2003) and (Davies and Foye, 1991)] for hopper flows but, these cannot be directly used to calculate the discharge rate in rectan gular hopper system examined in the present study. There have been a few studies on predicting th e cohesive discharge rate from hoppers, but no widely accepted correlati on exists for cohesive granular material. Segregation Segregation phenomena are widely present in ha ndling of granular ma terials, and authors such as Ottino and Khakhar (2000) and Herrmann and Luding (1998) have reviewed some of these phenomena and efforts to understand them. Se veral mechanisms for segregation have been proposed in the lite rature; de Silva et al (2000) summarized these mechanisms as follows: trajectory, air current, rolling, sieving (or sifting), impact, embedding, angle of repose, pushaway, displacement, percolation, fluidization, agglomeration, and concentrationdriven displacement. These causes of segregation have been subsequently condensed to a few key mechanisms [(Carson et al., 1986) and (Johanson, 1996)] by assuming that some occur infrequently such as embedding or pushaway and that others can be classified as a subset of another mechanism such as sifting and displaceme nt being special cases of percolation. These primary segregation mechanisms are summarized as: percolation, or sifting, segr egation whereby there is a net movement of smaller particles through a matrix of larger par ticles in the direction of grav itational acceleration, and that is enhanced by local strain, angle of repose segregation whereby material with a smaller angle of repose flows over material with a larger angle of repose during, for example, filling of a bin or silo, 19 PAGE 20 trajectory segregation whereby frictional effect s more effectively reduce the velocities of smaller or more angular particles in a thin layer of material moving down a chute, which results in differing trajectories upon discharge from the chute, and fluidization, or elutriation, se gregation whereby differences in drag forces on particles in a mixture will tend fluidize or entrain the fine particles causing them to segregate. Of these four mechanisms, it is assumed in the present work that percolation is the predominant means by which segregation will occur during hopper discharge for the following reasons. First, the granular material consists of relatively large, freeflow ing, spherical particles with uniform characteristics. Therefore, differe nces in angle of repose should be relatively small. Additionally, material discharging from a hopper generally does not experience projectile motion while within the hopper. Finally, it is assu med that intersti tial air effects are negligible due to the relatively large particle size. W ith these assumptions, the percolation mechanism should be the most prevalent in this system. Ketterhagen et al. (2008) did a thorough st udy of cohesionless segregation in a 3D wedge shaped hopper. This work extends that work for segregation of cohesive particles. Cohesion Granular materials are commonly used in the chemical, food and pharmaceutical industries. These materials are frequently produced by adding a liquid phase to mixtures of fine powders so that the liquid binds the smaller primary particles together into larger agglomerates. The liquid phase is then removed by convective or evaporative drying to yield an agglomerated product with improved properties that is suitable for commercial sale or further processing. Such agglomerated (or granulated) materials are t ypically denser and more freeflowing than a simple powder blend, tend to be less dusty, convey more consistently, and have improved volumetric filling and compre ssion characteristics. 20 PAGE 21 Granular m aterials have to be handled and conveyed in both the wet and dry state in a manner that ensures that they flow readily and do not segregate. For wet inprocess materials, this presents a big challenge as they are typically highly cohesi ve and the standard approaches for predicting flow and segregation performan ce of cohesionless powders are not directly applicable. Cohesion is the attractive force between particle s, and if cohesion is significant, substantial deviation from the freeflowing be havior of particulate systems is evident. Cohesion between particles can originate from several sources includ ing van der Waals forces, electrostatic forces, and liquid bridges (capillary forces). For part icles larger than several hundred microns, which are also free of electrostatic charge, the most important sour ce of cohesion is liquid bridges (Seville et al. 2000). This cohesion force is the focus of this study. Objectives This work focuses on simulating granular fl ow in a quasi3D, wedgeshaped hopper using the Discrete Element Method (DEM). The particle s are modeled as soft, inelastic, frictional, spheres. Figure 11 represents the outline of the computational work. The computational work can be divided into four parts: cohesionless discharge rate, cohe sive discharge rate, cohesionless segregation and cohesive segrega tion. Out of these, Ketterhegan et al (2008) did extensive work on the cohesionless segregation in a wedgeshape d hopper, which will not be discussed in this thesis. This work extends that work to cohesi ve segregation that will be discussed here. Numerous parameters that may affect the di scharge rate are studied. These parameters include hopper angle, hopper diamet er, particlehopper wall friction, particleparticle friction, fill height, particleparticle coeffici ent of restitution, and particle size. The predicted hopper discharge rates using DEM simulations are comp ared to the discharge rates obtained from empirical correlations available in the literatu re. In this way, the accuracy of the DEM 21 PAGE 22 sim ulations can be assessed and the sensitivity of the various parameters on the discharge rate can be explored. The pendular liquid bridge force is implemente d into the DEM code to model the cohesion between the particles. The stre ngth of the cohesive bond is ch aracterized by a Bond number. The effect of cohesion on discharge ra te is investigated. Specifically, four important parameters Bond number, liquid content, outlet width and hopper a ngle are studied to see how they affect the cohesive discharge rate. One of the goals of this study is to develop a co rrelation for predicting discharge rate for cohesi ve granular material. In a bidisperse system, the effect of cohesion on granular segregation is investigated. Bond number and liquid content are vari ed to study the effect on segregat ion. The goal is to find a set of parameters that minimizes segregation. Finally, three different charac teristics of hopper flow: discharge rate, angle of repose of material left in the hopper and size of clumps formed during cohesive discharge are investigated experimentally to validate the model and to better understand hopper discharge phenomenon and granular flows. 22 PAGE 23 23 Figure 11. Outline of computational work PAGE 24 CHAP TER 2 COMPUTATIONAL MODEL Introduction This chapter describes the computational model used in this work. First, a brief overview of DEM simulation is provided followed by the c ontact model used and th e description of the computational domain used. Then, the pendular liquid bridge force model used to include cohesion in this work described. Finally, the simulation parameters used are listed. DEM Simulations Discrete Element Method (DEM) simulation ha s emerged as one of the most important recent tools in probing granular flows. The method is extremely general in that Newtons laws of motion are used to determine the trajectories of individual particles. The time evolution of these trajectories then determines the global flow of the system. At the onset of a simulation, all particle prope rties are defined from a specified density and particle size distribution. The hoppe r geometry and dimensions are also defined. The particles are then inserted into the computational dom ain by defining a position and velocity for each. The simulation proceeds to a contact detection stag e where all particlepar ticle and pa rticlewall contacts are identified. For each contact, the deformation of the particles, which is modeled as an overlap, is calculated. Using the overlap va lues, forcedisplacement relations are used to calculate the forces acting on each particle. The total forces and moments acting on each particle are then summed and Newtons Second Law is used to calculate the translational and rotational accelerations. The time step is incremented and th e accelerations are integr ated over time using Eulers method to determine updated particle ve locities and positions. After measurements are taken, the process repeats until some ending condition, full discharge from the hopper for 24 PAGE 25 exam ple, is met. The key elements within the DEM algorithm are further detailed in the following sections. Contact Detection At every time step in a DEM model, all co ntacts including both particleparticle and particlewall interactions must be identified. Th is process is split into two phases and is called a nearneighbor search algorithm In the first phase, a list of all the particles in the neighborhood of each particle is created. Having created this list, in the second phase, the actual contacts and overlaps of each particle with all the other particles only in its neighborhood list is checked. This neighborhood list is then updated only after a certain number of time steps depending on its size. This algorithm helps save a lot of computati onal time compared to the brute force method in which each particle is checked for contact wi th every other particle at each time step. Contact Model There are two contact models commonly used wh ich are referred to as the hard or softparticle approaches. The hardparti cle approach assumes that particle s are rigid so that collisions are instantaneous and binary. As a result, hardp article models are generally best suited for dilute, collisional flows where thes e assumptions are good approximations. The soft particle model which allows for particle deformation and is modeled as an overlap of the spheres was used to calculate the forces between the particles in our system. The amount of overlap is a function of time step used. By using small enough time step, the overlap is generally maintained to 0.1.0% of the particle size. The softparticle model relies on a forcedisplacement (and/or forcedisplacement rate) relation to determine the interaction forces between particles. This approach proceeds via small time steps and is referred to as timedriven. Softparticle methods are generally utilized for de nse, enduringcontact flows as in our system. 25 PAGE 26 Accurate integration of the resu lting particle equations of m oti on dictates a small simulation time step and, hence, long computation times. For spherical particles, the overlap is calculated as 1 () 2jijiddxx (21) where d is the particle diameter and x is the position vector of the particle center for the respective particles i and j as shown in Figure 21. The contact forces between the particles ar e calculated using the model developed by Walton and Braun (1986).The normal force is calcula ted as: (22) where and are the loading and unloading sp ring constants, respectively, is the overlap between particles, and 0 is the overlap at which the unloading force is zero due to plastic deformation of the particles. The unit normal vector connecting the cen ter of two spherical particles is The effective coefficient of restitution for the contact, e is given by 0 for loading for unloading,L N Uk k n F nLkUk n .L Uk e k (23) The tangential force is calculated as: min,SN Tk FF s (24) where is the coefficient of sliding friction, is the tangen tial spring stiffness, and Tk is the tangential displacement (Cundall and Strack, 1979). The unit tangential vector is denoted by The f orce given in Eq. 23 is limited by the Coulomb criterion which specifies that the maximum s 26 PAGE 27 f rictional force is proportional to the normal force, where is the proportionality constant. The tangential displacement is determined by integrating the re lative velocity, x at the contact point: 0 ''t ttdt x s (25) where t0 is the time at which contact was initiated and t is the current time. The maximum value of the tangential displacement befo re gross sliding takes place is max/NFkT (26) Various contact models, includi ng the contact model used in the present work, have been recently compared (KruggelEmden et al., 2008). The results of this study indi cate that all of the tangential force models investig ated produce final collision prope rties that compare well with experimental results for most material combinat ions. In addition, the effects of the normal contact model on the results in th e tangential direction are small. A rolling friction, or rolling resistance model is also included in an effort to account for the nonspherical nature of real materials. This rolling resistance is modeled as an applied torque MR: Ri RN i MF (27) where R is the coefficient of rolling friction with units of length and is the angular velocity of particle i (Z hou et al. 2005). This moment acts in the dire ction opposite to the particle rotation thereby reducing the rotational velocity of a particle. i Once the co ntact forces for each contact are ca lculated, the body force due to gravitational acceleration FG is calculated as GiiFm g (28) 27 PAGE 28 where g is the gravitatio nal acceleration. Besides these forces, the other forces acting on a wet particle are the liquid bridge fore and viscous force which is disc ussed in the next section. Computational Domain The computational domain is a quasi3D plane flow hopper with periodic boundary conditions used on the front and back wall as shown in Fig.22. The hopper has rectangular crosssection with the lateral walls modeled as planar surfaces. Th e use of periodic boundary conditions greatly increases the computational effi ciency and, hence, only a thin slice of the hopper needs to be simulated to predict the behavior of a fully 3D system. The hopper is defined by the following ge ometric variables: the hopper width, W hopper outlet width, hopper depth, Zdepth, the hopper angle as measured from the vertical,0W and the initial fill height, H. Liquid Bridge Model When a quantity of liquid is introduced betw een particles, various regimes of liquid distribution are possible, depe nding upon the degree of saturatio n; these are the pendular, funicular, capillary and droplet regimes. For low liquid content, pendular liquid bridges are formed as shown in Fig. 23. Pendular liquid bridges are such that the interaction betwee n the particles is always binary. As the liquid content increases, the liquid bridges are not inde pendent from each other and the pendular liquid bridges begin to merge as the system transitions to the funicular and capillary regimes. At very high liquid content, the droplet regime is reached and the particles are essentially immersed in the liquid. For the present study, we only consider the pendular regime. The stable, smooth liquid bridge formed between two spherical bodi es results in an attractive force between them. Th is force, in the case where grav itational distortion of the liquid 28 PAGE 29 bridge is negligible, contains two components: (1) th e axial com ponent of the surface tension acting on the threephase contact li ne, and (2) the hydrostatic force due to the pressure deficiency in the bulk of the liquid. The theoretical analysis of such a liquid bridge force has been the focus of previous studies (see, for example, Mehrotra and Sastry (1980), Lian et al (1993) and Willett et al (2000)). Lian et al (1993) derive the liquid bridge force as a function of the surf ace tension, halffilling angle, contact angle, and particle radius. However, to implement this force model into a DEM code in a computationally efficient manner, an explicit relationship between the force, liquid bridge volume, and contact angle and particle separa tion is needed. The regression expressions generated by Mikami et al (1998) and Soulie et al (2006) provide just this sort of relation. Mikami et al (1998) developed a dimensionless cohesive force correlation for a pendular bridge spanning two equalsized spheres as a function of the liquid volume and particle separation distance. Soulie et al (2006) extended this expression for two unequalsized spheres as: 12 2exp()lbh FRRAB R C (29) 0.53 3 21.1() V A R (210) 2 3 22(0.148ln()0.96)0.0082ln()0.48 V B RR 3V (211) 3 20.0018ln()0.078 V C R (212) where is the surface tension of the liquid, 1 R and 2 R are the radii of the two spheres, is the contact angle of the liquid with the particle, h is the separation distance, and V is the liquid 29 PAGE 30 volum e associated with the liquid bridge. The constants A B and C are dimensionless regression parameters. Another important parameter for the liquid bri dge force between particles is the rupture distance hrupture. The rupture distance is the maximum se paration distance between the particles for which the liquid bridge is stable. Experime nts by Mason and Clark (1965) show that the rupture distance for any meridianal profile of the liquid bridge vari es linearly with contact angle and V1/3. Additional theoretical cons iderations as given by Lian et al (1993) have yielded the following specific relationship between the rupture distance and the liquid bridge volume and the contact angle: 1 3(10.5)ruptureh V (213) In addition to the liquid bridge force, partic les are subject to a vi scous, resistance force Fvf resulting from the squeezing out and pulling in of liquid between th e two closely spaced particle surfaces. Assuming rigid spheres, the viscous sq ueeze film force in the normal direction is: *6vf n R FRv h (214) 12111 R RR (215) where is the viscosity of the liquid, vn is relative speed between the particles, and R* is the effective particle radius (Lian et al.,1998). A minimum value for h, corresponding to the combined asperity height of the tw o surfaces, in this case equaling 1/100th of the smaller particle diameter, is employed to keep the viscous force b ounded. It has been observed in this work and in the work by McCarthy (2003) that segregatio n and mixing results are insensitive to a two order change in the magnitude of the chosen aspe rity value. A similar tangential viscous force 30 PAGE 31 m ay also be derived, but previous studies have shown that this force ha s a negligible effect on particle motion (Hsiau and Yang, 2003). An important parameter for characterizing th e relative cohesiveness of a system is the Bond number, Bo which is defined as the ratio of the maximum cohesive force acting on a particle, Fc, to the gravitational force acting on the particle, W : 2 32 4 2 3CF R Bo WR Rg 3 g (216) where is the surface tension of the liquid, R is the radius of the particle, is the density of the particle, and g is gravitational acceleration (Nase et al. 2001). There are two Bond numbers associated with a bidisperse system: Bol which corresponds to the large particles, and Bos which corresponds to the smaller particle s. As will become evident later, the Bond number associated with the smaller particles, i.e. Bos, is used to characterize cohesi on throughout this work. As the Bond number of the system increases, the cohesi veness of the system in creases. For example, Bo = 0.25 represents a maximum cohesive force on a smaller particle to be 25 % of the smaller particles weight. Assumptions In implementing the liquid bridge force model into the DEM code, the following assumptions regarding liquid distribution, bond formation, and bond rupture are used. 1. The total liquid content of the system is distributed so that every particle has a constant thickness of liquid layer around it. Hence, in a bidisperse sy stem the ratio of the total amount of liquid on a larger particle to the liquid on a smaller particle is (dL/ds) 2. 2. When two particles come into contact, a liquid bridge is formed between them. 3. When particles are in an enduring contact, a Hookean spring force, a Coulombic friction force, and an attractive cohesive liquid bridge force act on the particles. The magnitude of 31 PAGE 32 the cohesive force rem ains constant during th e contact as shown sc hematically by the arm AB of Fig. 24. 4. When particles are separating afte r a contact, only gravity, a liquid bridge force and a viscous force in the normal direction act. The liquid br idge force for this situation is represented schematically by arm BC in Fig. 24. 5. The amount of liquid from each particle that take s part in bond formation is determined using the scheme proposed by Shi and McCarthy (2008). The liquid on the surface of each particle that is within the area of the spherical cap neighbori ng the contact spot will contribute to the liquid bridge. In other words, the bridge harve sts liquid from the particle surfaces near the contact spot. With two part icles contacting each other, the spherical cap is bounded by tangential lines from the cente r of one particle to the su rface of the other particle. 6. If the separation distance between particles ex ceeds the critical rupture distance for a given liquid volume, the liquid bridge ruptures and the liquid returns to each of the particles. It is assumed that since the contact angle of the liquid with both particles is the same, the liquid transfer between the particles is negligible. Hence, each particle car ries the same initial amount of liquid during the cour se of entire simulation. 7. Liquid bridges slip over the particle surface; that is, there is no attractive force in the tangential direction due to liquid bridges (t here is only a normal attractive force). 8. Liquid bridges that exist between a particle and the hopper wall are treated in the same manner as liquid bridges between particles. Al so, all of the liquid in a particlewall bond comes from the particle. The wall is assumed to be dry. Simulation Parameters The computational work in this thesis can be divided into three parts: the study of cohesionless discharge rate, th e study of cohesive discharge ra te and the study of cohesive segregation. A set of baseline parameters was established for each of these sets of DEM simulations. All the computational work employ these parame ters for each set of DEM Simulations unless otherwise noted. In order to study the effect of one individual parameter, its value is modified from its baseline value. Table 21 lists the ba seline parameters that were used for the DEM simulations to study cohe sionless discharge rate. 32 PAGE 33 The particle s modeled in the simulations are treated as glass spheres and hence a particle density of 2.5 g/cm3 is assumed. The value of was based on dense phase shear sim ulations which showed that the particle stress is insensitive to provided it is sufficiently large (Ketterhagen et al. 2005). The ratio determ ines the coefficient of restitution. For glass spheres, the coefficient of restitution of 0.94 is used for particleparticle collisions, whereas a value of 0.90 is used for particlewall collis ions. It is also important to note that determines the period of normal oscillation and determines the period of tangential oscillation. For uniform spheres, the periods of tangential and normal oscillation are equal when LkLk /LUkkTkLk /2/7TLkk (Cundall and Strack, 1979). A smal l rolling resistance value is commonly used to take into account surface asperities and/or a slightly nonspherical shape, which act to inhibit rolling in real materials. However, the discharge rate results are insensitive to this parameter. In our simulations a value of 0.01cm is used. The use of periodic boundary conditions reduces the size of the computational domain, and thus the number of particles and computational time. It is desired to use the smallest computational domain that still produces the same results as a system without periodic boundaries. In the present work, the discharge rate is independent of depthzd for values greater than 2.1The hopper width is chosen to be suf ficiently large so that it does not affect the discharge rate. Both Ketchum (1929) and Brown and Ri chards (1960) suggested that as long as W / W0 > 2.5, the discharge rate will be independent of hopper width. Hence, the baseline condition for dimensionless hopper width was chosen to be W / W0 = 5. The reason for choosing the baseline 33 PAGE 34 conditions for the fill heighttohopper width ratio and the fricti on coefficient will be discussed in the Resu lts section. Table 22 lists the range of various hopper and particle properties that were used in developing the various experiment al hopper correlations, along with those used in the current DEM simulation work. The baseli ne conditions, shown in Table 21, lie within this range. In order to study the cohesive discharge rate identical hopper paramete rs are used as in the cohesionless discharge rate. Although, the partic les modeled in the simulations are treated as cast iron spheres and hence a particle density of 7.85 g/cm3 is assumed. The liquid for the base case is assumed to be water and the base case liq uid content is taken as 2% by volume. This system has a Bond number of 1. The Bond number is varied by changing the surface tension of the liquid, but the viscosity and contact angles of the liquid with different surfaces is assumed to be identical to that of pure wa ter. Table 23 lists the baseline parameters for the cohesive discharge rate. Cohesionless segregation in a hopper was studied by Ketterhagen et al. (2008) in detail. To study the cohesive segregation id entical hopper parameters were used. The particle diameter ratio, D hopper angle, and initial fines mass fraction, xf, were chosen such that discharge segregation for the case of cohesionless par ticles was maximized for these parameters (Ketterhagen et al. 2008). Table 24 lists the baseline para meters for the cohesive segregation. 34 PAGE 35 j iLkULkk ir Oi x j x x y z Sk n s Figure 21. A schematic of normal and tangential contact force models 35 PAGE 36 zdepthW H Wo z y x Figure 22 Schematic of the quasi3D computational domain that models a thin slice of a large, rectangular crosssectioned hopper. 36 PAGE 37 Figure 23. Representation of the liquid brid ge formed between particles of unequal size Figure 24. Illustration of the liquid bridge force as a function of the separation distance between two particles. 37 PAGE 38 Table 21. L ist of baseline parameters for study of cohesionless disharge rate. Particle Properties Total number, N 6790 Mean diameter 0.224 cm Density 2.5 g/cm3 Interaction Parameters Spring stiffness particleparticle loading kL 250,000 g/s2 unloading kU 283,000 g/s2 tangential kT/ kL 2/7 particlewall loading kL,w 250,000 g/s2 unloading kU,w 308,600 g/s2 tangential kT,w/ kL,w 2/7 Friction coefficient particleparticle sliding 0.20 rolling R 0.01 cm particlewall sliding w 0.20 rolling R,w 0.01 cm Dimensionless hopper geometry Aspect ratio, A = H / W 0.74 Width, W / Wo 5.0 Outlet width, Wo/ d 11.2 Depth, Zdepth 0.56 cm (Zdepth/d = 2.5 for base case) Halfangle from vertical, 90 deg 38 PAGE 39 Table 22. L ist of hopper and particle properties used in experimental correlations Particle type Particle diameter ( d) Hopper outlet width ( Wo/d ) Aspect ratio ( H/W ) Hopper angle ( ) Beverloo et al. (1961) Linseed, Spinach, Watercress, rapeseed, Kale, Swede, Turnip Seeds Sand Wheat Crystallized sugar Polystyrene globules Green peas Soybeans Lupin 0.0093 cm (Sand) 0.7 cm (Green peas) 4.2 (Sand) 15.6 (Spinach) 0.667 6 90 Brown and Richards (1965) Beads Sand Tapioca Pentremawr Coal 0.0015 cm (Coal) 0.168 cm (Tapioca) 2.8 (Tapioca) 21.5 (Sand ) Rose and Tanaka (1959) Silica sand Ground glass Granulated sugar Steel balls 0.0112 cm (Sand) 0.318 cm (Steel balls) 5 (Steel balls) 104.82 (Sand) 15 90 Humby et al. (1998) ABS Acrylic Radish Turnip 0.06 cm (Acrylic) 0.325 cm (ABS) 7.4 (ABS) 39.6 (Acrylic) 4 90 Current DEM simulation work Glass spheres 0.116 cm 0.224 cm 6.7 24.6 0.74 2.96 25 90 39 PAGE 40 Table 23. L ist of baseline parameters for study of cohesive disharge rate. Particle Properties Total number, n 6790 Mean diameter, d 0.235 cm Density 7.85 g/cm3 Interaction Parameters Spring stiffness particleparticle loading kL 250,000 g/s2 unloading kU 283,000 g/s2 tangential kT/ kL 2/7 particlewall loading kL,w 250,000 g/s2 unloading kU,w 308,600 g/s2 tangential kT,w/ kL,w 2/7 Friction coefficient particleparticle sliding 0.20 rolling R 0.01 cm particlewall sliding w 0.20 rolling R,w 0.01 cm Cohesion viscosity of liquid, 0.001 g/cm/s contact angle, liquidparticle, 4 deg contact angle, liquidwall, 70 deg liquid content, L 2 % (by volume) surface tension, 72 dyn/cm Bond number Bo 1 Dimensionless hopper geometry Aspect ratio, A = H / W 0.74 Width, W / Wo 5.0 Outlet width, Wo/ d 11.2 Depth, Zdepth/ d 2.5 Halfangle from vertical, 90 deg 40 PAGE 41 Table 24. L ist of baseline parameters for study of cohesive segregation. Particle Properties Total number, n 33470 Mean diameter ( dL and dS) 0.224 cm and 0.0521 cm Density 2.5 g/cm3 Fines mass fraction, xf 5% Interaction Parameters Spring stiffness particleparticle loading kL 250,000 g/s2 unloading kU 283,000 g/s2 tangential kT/ kL 2/7 particlewall loading kL,w 250,000 g/s2 unloading kU,w 308,600 g/s2 tangential kT,w/ kL,w 2/7 Friction coefficient particleparticle sliding 0.20 rolling R 0.01 cm particlewall sliding w 0.20 rolling R,w 0.01 cm Cohesion viscosity of liquid, 0.001 g/cm/s contact angle, liquidparticle, 4 deg contact angle, liquidwall, 70 deg Bond number (smaller particles), Bo 1 liquid content, L 2 % (by volume) Dimensionless hopper geometry Aspect ratio, A = H / W 0.74 Width, W / Wo 5.0 Outlet width, Wo/ dL 11.2 Depth, Zdepth/ dL 2.5 Halfangle from vertical, 90 deg 41 PAGE 42 CHAP TER 3 EXPERIMENTAL Introduction This chapter discusses the details of the experi mental component of this thesis work. First, the motivation behind the experiments are discusse d .Then, the granular media and hoppers used in the experiments are described. Then, the pr ocedure for carrying out the experiments and the method to collect the data is presented in the next section. Motivation Most of the work described in this thesis is computational in nature, but a few basic experiments were carried out in hoppers of dimens ions similar to that of computational work to validate the computational model. Also, compari ng the visualization from the simulation and the results from experiments gives fu rther insights into the hopper di scharge phenomenon. Also, it is possible to get a better estimate of some of the parameters used in the simulations with the help of these experiments. For example, particleparti cle friction coefficient for the cast iron shots can be estimated. Materials and Experimental setup Two types of granular media are used in the experiments. Glass beads of average size 0.22 cm and cast iron shots of average size 0.235 cm. The granular media is relatively spherical and free flowing. In order to introdu ce cohesion into the system two t ypes of liquids are used with different surface tension. Water with a surface tens ion value of 73 dyn/cm and silicone oil with a lower surface tension of 20.8 dyn/ cm are used. The resulting four cohesive systems that are investigated experimentally are listed in Table 31 42 PAGE 43 In addition to the cohesive cases listed in Ta ble 31, cohesionless ca st iron shots and glass beads are als o used in the experiments. Cast ir on shots with silicone oil is the least cohesive system whereas glass beads with water is the mo st cohesive system. Cast iron shots with water and glass beads with silicone oil have iden tical bond number sugges ting they should show similar cohesive properties. The hoppers are shown in Figure 31 and are made of clear acrylic for visualization of the flow. The hoppers have dime nsions of 12.5cmx12.5cmx100cm. There is a 12.5cmx 2.5 cm rectangular orifice at the bottom equipped with a gate as shown in Fig. 31C. The hoppers have halfangles of 90 (flatbottomed) and 55 with respect to vertical. The Ohaus Champsquare balance was used with Collectv6.1 software to get the weight of the material remaining in the hopper. Procedure This section outlines the experimental pro cedure and the method to collect data. Three different characteristics of flow were investigated: discharge rate, angle of repose of the material left in the hopper and the clump si ze in the cohesive discharge. For calculating the discharge rate, the balanc e is connected to the computer using the Ohaus proprietary cable, and Collectv6.1 software is used to get the weight readings on the balance. At the start of each discharge rate experiment, the hopper is cleaned and placed on the Ohaus balance. For the case of cohesion, the appropriate liquid is uniformly mixed with the granular media of interest in a separate contai ner. These wet particles are then poured into the hopper. The collectv6.1 software records the wei ght readings on the balance every 0.4 seconds directly to an Excel sheet. Th e hopper slide gate is opened, and the particles coming out of the hopper are collected in a container that is not pl aced on the balance. At any point of time the 43 PAGE 44 balance read s the amount of material remaini ng in the hopper. By plotting the data of mass remaining in the hopper vs. time and calculating the sl ope of this plot the discharge rate can be calculated. A set of three experimental runs ar e made for each system. The sources of error in calculating the discharge rate are due to the nonu niformity of liquid distribution in the system, the sensitivity or the stabilization time of th e Ohaus balance in recording the data and in determining the region of constant discharge rate following the initial transient just after opening the hopper gate. For calculating the angle of re pose of the material left in the hopper, a snapshot of the cross section of the hopper is take n after the discharge is complete as shown in Fig. 32A for cohesionless glass beads. A set of two snapshots fr om each system is analyzed to get the value of average angle of repose. From this snapshot the angle of repose can be calculated. To calculate the angle of repose from the simulations a similar snapshot is taken from the visualization of the flow and is shown in Fig. 32B. The only source of error is in determini ng the free surface of the particles in highly cohesive system. For calculating the size of the clumps, a high sp eed camera at 400 fps is used to capture the particles coming out of the hopper. Individual frames are extracted from the video of the flow. In general, a distribution of clump sizes is obser ved for any given Bond number. The largest sized individual clumps which are generally formed n ear the end of discharge are visually identified and isolated. These raw images are then analyzed for the number of pixels within the largest clump. By knowing the number of pixels occupied by a single partic le, the total number of particles in the clump can be estimated. Figure 33 shows one of such clumps formed in the cast iron shots with water system. The number of pi xels occupied by the central large clump is 1498. One of the particles takes approximately 100 pixels, suggesting that the clump is made of 44 PAGE 45 roughly 15 particles. A set of three im ages are analyzed in this way to give the average number of particles in the largest clumps formed at any given Bond number. To calculate the size of clumps from the simulation, the visualization of the flow can be directly used in a similar manner to calculate the number of particles in a clump. A Figure 31. Schematic of the experimental hoppe r. A) Assembly B) Crosssection of the 55 hopper C) Crosssection of the 90 hopper D) Gate and orifice of the hopper 45 PAGE 46 B C Figure 31 continued. 46 PAGE 47 D 47 PAGE 48 A B Figure 32. Snapshots for calculating angle of repose. A) Experiments B) Simulation 48 PAGE 49 Figure 33. Example of a snapshot from hi gh speed camera for calculation of clump size Table 31. Types of experimental cohesive system System Bond Number Cast iron shots with silicone oil 0.3 Cast iron shots with water 1 Glass beads with silicone oil 1 Glass beads with water 3.5 49 PAGE 50 CHAP TER 4 RESULTS AND DISCUSSION Introduction This chapter describes the results obtained fr om the computational and experimental work. The first part of this chapter describes the computational results for discharge rate of cohesionless particles. The second part describes the computational results for discharge rate of cohesive particles. The third part describes the computational results for segregation of cohesive bidisperse particles. Finally, it presents the e xperimental results and compares them to the simulation results. Discharge Rate for Cohesionless Particles This section examines the effects of various pa rameters that may affect the discharge rate from a hopper. The results obtained from DEM simulations are compared to the various experimental correlations that we re discussed in Section 2.1. It should be noted that a value of 0.59 is used for the poured solid fraction of unshaken glass spheres (Cumberland and Crawford, 1987) when calculating the mass discharge ra te for the experimental correlations. Fill Height and Friction As has been discussed previously, the discha rge rate of coarse granular materials is generally assumed to be independent of the fill height. Fig. 41A shows the discharge profile when the friction coefficient is 0.05 (the particle particle and particlewal l friction coefficients are the same) for various values of fill heights. Fig. 41B shows a similar plot for a larger friction coefficient of 0.84. For smaller friction coefficients, the discharg e rate increases with increasing fill height, which is evident from the increasing slopes of the discharge profile (Fig. 41A). For the case of large friction coefficients, all the four discharge profile for the various values of fill height have 50 PAGE 51 the sam e slope indicating an independence of the discharge rate on fill height (Fig. 41B). These observations point to fluidlike behavior of bulk solids at very small friction coefficients. Based on this observation, a parametric study wa s performed to determine at which values of the friction coefficients and w the discharge rate becomes independent of fill height. Figs. 4a 4f show the results obtained from this parametric study. In these figures, the instantaneous discharge rate is plotted against the instantane ous fill height. Figs. 4a4f indicate that the discharge rate becomes independent of the fill height at a critical value of = w = 0.2. For friction coefficients less than 0.2, Figs. 42A 42C show that the mass discharge rate increases with increasing fill height. This behavior is similar to the behavior of a fluid. However, for friction coefficients greater than 0.2, the mass di scharge rate is independent of fill height, as shown in Figs. 42D 42F, which is consistent w ith Janssens theory. Hence, a critical value of 0.2 was chosen for the baseline condition for the friction coefficients. For cases in which the friction coefficients are varied fr om the baseline value, only fricti on values larger than 0.2 were used. Hence, all of the discharge rate results reported henceforth are independent of fill height. It is postulated that most of the naturally occurring materials like corn which were used in the published experimental studies like those of Janssen (1895) had a pa rticleparticle friction coefficient of at least 0.2 or larger. This postulate is supported by the work of Suiker and Fleck (2004) who observed a saturation of macroscopic friction angle of granular materials when the particleparticle friction coe fficient is greater than 0.2 For the case of low friction solids, the discharge rate will be a function of fill height. A hypothetical frictionless solid will behave like a fluid showing a fill height dependency of H0.5. 51 PAGE 52 ParticleParticle and ParticleWall Friction As has already been shown in the previous s ection, friction plays an important role in determining the hopper discharge rate In order to determine whic h type of particle friction particleparticle or particlewall (if any) domin ates the discharge rate, each friction coefficient was varied independently as summarized in Table 41. By comparing the results for the four cases, it is clear that particleparticle friction plays a much more important role than particlewall fr iction on the hopper discha rge rate. For the case of large particleparticle friction, the discharge ra te is approximately the same as the value for the large particleparticle and particlewall friction case. The large particlewall friction case has a discharge rate comparable to the base case (Fig. 43). These results suggest that particleparticle friction is the primary controlling factor in dete rmining the discharge rate with the particlewall friction coefficient having little effect. This finding can be explained in terms of the flow regime. In a broad hopper (in this case, W / W0 = 5), particlewall friction does not affect the flow regime near the outlet since ther e is typically a stagnant region located between the flowing region and the walls. Also, since most of th e interactions in a fl owing hopper system are particleparticle interactions, particleparticle fric tion is expected to play a more significant role in determining the flow regime and hence the discharge rate. Coefficient of Restitution The coefficient of restitution is a parameter which has not been studied in any of the experimental investigations to date. Fig. 44 plots the discharg e profile for two DEM simulations with significantly different coefficients of restitu tion (note that the partic leparticle and particlewall coefficients of restitution are identical). As is evident from the plot, changing the coefficient of restitution has negl igible effect on the discharge ra te. In another simulation study, 52 PAGE 53 Ristow (1997) observed that cha nging the coefficient of restitution from 0.3 to 0.7 showed a change in the discha rge rate of 1.2 %. The flow through a hopper is dominated by long dur ation, frictional, sl iding contacts. In such flows, the coefficient of restitution will ha ve a negligible influence on the flow dynamics. Hence, the friction coefficient is expected to ha ve a much more significant role in the system dynamics than the coefficient of restitution. Hopper Width Many experimental investigations have shown that hopper width has no effect on the discharge rate (Nedderman et al ., 1982). Fig. 45 plots the disc harge rate obtained through DEM simulations for two hoppers one twice the width of the other. Identic al discharge rates are observed, with the doublewidth hopper discharging over a longer pe riod of time sin ce it initially contains twice the amount of material. Thes e results are consistent with experimental observations. The reason for the behavior is likely due to the fact that as long as the hopper width is sufficiently large, the walls do not influence the flow re gime near the hopper outlet and, thus, have no effect on the mass discharge rate. Hopper Angle Hopper angle (from the vertical) clearly affects th e flow behavior in the hopper and its discharge rate. Fig. 46 plots the mass discharge rate as a function of hopper angle in addition to the predictions from the empirical co rrelations by Brown and Richards (1965) and Rose and Tanaka (1959). There is a good agreement between the DEM discharge rate and the predictions from the correlations. The discharge rate shows little variation as the hopper angle decreases from 90 to approximately 55 Beyond this critical angle, the discha rge rate increases rapidly as the hopper angle decreases. This critical angle is the fl ow angle of the stagnantflowing boundary angle for the material inside the hopper. Th ere is little or no change in the funnel flow pattern within the 53 PAGE 54 hopper until the hopper angle is less than this critical angle. As soon as the hopper angle decreases below this angle, the m aterial starts flowing freely inside the hopper in a mass flow manner and the discharge rate increases. The Ro se and Tanaka (1959) correlation seems to fit the DEM discharge rate data better for angles below the critical angle, whereas Brown and Richards (1965) correlation works better for an gles larger than the critical angle. Outlet Width Outlet width is reported as being the mo st significant parameter for controlling the discharge rate from a hopper. Dimensional anal ysis suggests that the mass discharge rate is a function of outlet width raised to the 5/2 power for circular or ifices. The modified Beverloo equation for rectangular hoppers dire ctly predicts the discharge rate as a function of outlet width (corrected for dead space: Z=kd ) raised to the power 3/2 as show n in Eq.13. Fig. 47 compares the mass discharge rates obtained from the DEM simulations to predictions from the Beverloo correlation (Eq. 13). The mass discharge rate pe r unit length is raised to the 2/3 power and plotted against the outlet width. The straight line is the trend us ing the Beverloo correlation. As can be observed, the DEM simulations match the Beverloo trend extremely well. The intercept of the straight line on the abscissa gives the Z = kd term in the Beverloo correlation. Knowing the particle diameter d, k can be calculated. A value of 1.3 is obtained from the simulations, which is close to the value of 1.5 0.1 predicted by experiments for spherical particles (Nedderman and Laohakul, 1980). Particle Diameter Particle diameter is expected to have only a minimal influence on the discharge rate as suggested by the Beverloo correlation. Fig. 48 shows that the discharge rate varies only 3% with a 30% decrease in particle diameter. For a particle with a smaller diameter, the Z = kd term will be smaller indicating that the effective flow ing region exit area will be larger and a larger 54 PAGE 55 flow rate will result. The Beverloo correlation pr edicts a change of approxim ately 5% in this case Particle Size Distribution Humby et al (1998) extended the Beverloo correlation to account for bidisperse particulate systems (Eq. 16 and 17). The correlation indi cates that the discharge rate increases with increasing fines mass fraction. The rate at whic h the discharge rate increases with fines mass fraction decreases as the fines mass fraction increases, as can be obser ved in Fig. 49. This trend is the result of the increasing flowing density of th e particulate material as the fines fill the voids between the coarser particles (Dias et al., 2004). This filling of voids between coarser particles with fine particles increases the flowing bulk density until all the voids are filled with fine particles. After this point there is a li ttle effect of adding fines to the system. Bidisperse particle DEM simulations were also performed for different values of fines mass fraction with coarse particles having the ba seline diameter of 0.224 cm and fines having a diameter of 0.116 cm. The DEM discharge rates, shown in Fig. 49, compare well with the Humby correlation; however, the DEM values are on an average 8% smaller than the values determined from the Humby correlation. Though this difference is small, it is consistent. The cause for the discrepancy remains unclear. Discharge Rate for Cohesive Particles This section examines the effect of four im portant parameters on the discharge rate of cohesive particles from rectangular hoppers. These parameters are the outlet width, Bond number, liquid content and hopper angle. Previous simulation studies for cohesionless particles indicate that particle diameter, hopper width and coefficient of restitution have minimal effect on discharge rate (Anand et al., 2008), and hence these parameters were not studied here. Based on the simulation results, a Beverlootype correlatio n may be developed for predicting the discharge 55 PAGE 56 rate from the flatbottomed hopper by modifying th e Beverloo constant for cohesive particles as a function of Bond number. Effect of Outlet Wi dth and Bond Number As mentioned before, the Bond number character izes the cohesiveness of particles. A value of Bo = 0 represents cohesionless particles, whereas Bo = 1 represents a system in which the maximum cohesive force acting on a pa rticle is equal to its weight. In order to understand the difference between cohesive and cohesionless systems, several simulations in hoppers of different dimensionless outlet widths ( Wo/d = 6.4, 10.6, 15, 19.2 and 23.4) were carried out for varying Bond numbers (0 5). A dimensionless discharge rate for cohesive particles was defined by comparing it to th e discharge rate of cohesionless particles in a hopper of the same dimensions as follows: Discharge rate for cohesive system Dimensionless discharge rate = Discharge rate for cohesionless system (41) Fig. 410 plots the dimensionless discharg e rate as a function of Bond number and dimensionless outlet width ( W0/d ). For a fixed outlet width, the discharge rate decreases as Bond number increases. It is evident that the deviation from the cohe sionless case is dependent on two parameters: Bond number and dimensionless outlet width. These two parameters seem to be important in determining the regime of cohesive flow, and hence the devi ation in the discharge rate of cohesive particles when co mpared to the cohesionless case. The Beverloo equation for rectangular hoppers of exit width and length applied to cohesionless particles, predicts that the m ass discharge rate is proportional to the outlet width (corrected for dead space Z, Z=kd where k is the Beverloo constant and d the particle diameter) raised to the power 3/2 as: OWOL 56 PAGE 57 13 21.03()()flatbottom flowO OWg L k d W2k d (42) where f low is the flowing bulk density and g is the acceleration due to gravity (Myers and Sellers, 1971). Beverloo took dead space into account for the effects of a layer of particles just near the edge of the orifice wh ere there is little or no flow. Fig. 411 plots the mass discharge rate per un it length raised to the 2/3 power against the outlet width for various Bond number s. For the cohesionless case, it has been observed that the DEM predictions for the discharge rate match extremely well with the predictions from the Beverloo correlation yielding line ar behavior (Anand et al., 2008). Similar plots for cohesive particles (Bo = 1, Bo = 2 and Bo = 3) reveal th at the basic discharge ra te trend given in the Beverloo correlation, where the ma ss discharge rate varies linearly with the (corrected outlet width)3/2 is also valid in the case of cohesive particles. The intercept of the straight line on the abscissa of Fig. 411 gives the Z = kd term in the Beverloo correlation. Knowing the particle diameter d, k can be calculated. Fig. 412 plots the calculated value of k as function of Bond number. By regression analysis, a simple correlation of this Beverloo constant may be generated for different values of Bond number: 0.39 = 1.9 B oke. (43) This expression can be used to modify the Beverloo correlation, originally developed for cohesionless particles, for predicting the discharge rate of cohesive particles. This modification is physically justified by the fo rmation of visible clumps of cohesive particles during the discharge from the hopper. For th e cohesionless case shown in Fig. 413A, there is little dead space at the edges of the outlet, while in the cohesive case shown in Fig. 413B, the clumps tend to decrease the effective ou tlet width of the hopper by incr easing the dead space, Z. 57 PAGE 58 Fig. 414 compares predictions for the discharge rate obtained from the DEM simulations with the predictions from a Beverloo correlation modified for cohesion. 13 0.390.39 221.03[(1.9)][(1.9)]Bo Bo flatbottom flowO OWg Le d Wed (44) The average deviation b etween the two predictions for discharge rate for Wo/d = 6.4 and 10.6 are approximately 8% and 10%, respectively. The cohesion correlation does not capture the degree of curvature observed in the DEM data. One source of error in the DEM data arises from estimating the discharge rate when the flow is hi ghly cohesive and intermittent, especially near the condition of no flow. Another important source of error is the assumption that f low is a constant in the application of the Beverloo corre lation with varying cohe sion. A modification to f low taking cohesion into account will lik ely yield more accurate results. Although it is possible to generate alternate expressions for the discharge rate which fit the simulation data better, the simplicty of modifying only one parameter in the well established Beverloo corrrelation is compelling. Effect of Liquid Content The liquid bridge force is a weak functi on of the liquid content; however, the liquid content (refer to Eq. 213) direct ly affects the range of particle separation distances over which the cohesive force acts. Increasing the liquid content increases the rupture distance and, hence the distance over which cohesion acts in the system. It is important to note that at extremely small liquid contents, some of the assumptions made in the liquid bridge model will become in accurate. Specifically, the distribution of the liquid on a particles surface is unlikely to form a uniformly thin film that will be available for the formation of a liquid bridge for an arbitrary contact point on the partic les surface. Since the lower limit at which the distributed thin film assumption becomes inaccurate is unknown, an 58 PAGE 59 arbitra rily chosen lower value of L = 1% (by volume) was used. A maximum value of L = 23% (by volume) was chosen since this value corr esponds to the upper limit for pendular liquid bridges for a monodisperse system of spherical particles (Yang et al., 2003). These liquid contents are also within the ra nge of validity of the Mikami et al. (1998) regression expressions. Fig. 415A and Fig. 415B plot the mass disc harged from the hopper as a function of time for Bo = 0.3 and Bo = 2.0 respectively. The slope of this curve gives the mass discharge rate. At Bo = 0.3, liquid content plays a negl igible role in determining the discharge rate as evident from the identical slopes of the discharg e profile for different values of the liquid content. There is only a 2% variation in the discharg e rate over the range of liquid contents at this small value of Bond number. For Bo = 2, the effect of liquid content on the discharge rate becomes more important. At this value of Bond number, there is a decrease in the flow rate with increasing liquid content as observed by th e decreasing slope of the disc harge profile. There is a 14% variation in the discharge rate at this large value of Bond number over the range of liquid contents. Even at the larger Bond number the effect of liquid content on discharge rate is much less than the effect of the Bond number itself Changing the Bond number by a factor of approximately seven, i.e. from 0.3 to 2, causes a 33% decrease in the discharge rate, whereas for Bo = 2, increasing the liquid content by a factor of eight, i.e. from 1% to 8%, causes only a 6% decrease in the discharge rate. An important conclusion from this liquid bridge force model is that Bond number, which characterizes the strength of cohesive bonds, plays a dominating role in determining the discharge rate when compared to liquid conten t, which characterizes the range of the liquid bridge force. 59 PAGE 60 Effect of Hopper Angle Hopper angle (from the vertical ) strongly affects the flow be havior in the hopper and its discharge rate. Fig. 416 plot s the mass discharge rate as a f unction of hopper angle for different Bond numbers. Previous work (Anand et al., 2008) compared the DEM resu lts to the predictions from empirical correlations by Brown and Richar ds (1965) and Rose and Tanaka (1959). Good agreement was obtained between the DEM disc harge rate and the predictions from the correlations. For the cohesionless case (Bo = 0), th e discharge rate shows little variation as the hopper angle decreases from 90 to approximately 55 Beyond this critical angle, the discharge rate increases rapidly as the hopper angle decreases. This critical angle is the flow angle or the stagnantflowing boundary for the material inside th e hopper. There is little or no change in the funnel flow pattern within the hoppe r until the hopper angle is less than this critical angle. As soon as the hopper angle decreases below this angle, the material starts flowing freely inside the hopper in a mass flow manner and th e discharge rate increases. The data for cohesive systems follows the same trend, but the critical angl e decreases as the cohesion increases. This observation is consistent with the increase in th e angle of repose of cohesive materials when compared to cohesionless materials. Segregation of Cohesive Particles This section examines the role that various cohesion parameters play on the extent of discharge segregation from the hopper. The cohesive model can be characterized in terms of two important parameters. The Bond number determin es the strength of cohesive force acting between the particles, and the rupture distance de termines the range over which the cohesive force acts. Since the rupture distance is directly related to the volume of the liquid bridge and the contact angle as given by Eq. 211, the amount of liquid in a specific fluidparticle system directly affects the range over which the cohesive force acts. 60 PAGE 61 Effect of Bond Number on Segregation As the Bond number increases, the cohesiveness of the particles increases. A value of Bo = 0 represents cohesionle ss particles, whereas Bo = 1 represents a system in which the maximum cohesive force acting on the particle is equal to the particles weight. Fig. 417A plots the mass fraction of fines in a given hopper discharge sample,i x normalized by the mass fraction of fines initially in the system, f x as a function of the fractional mass discharged from the hopper, TotalMM for Bo < 1 and a liquid content of L = 2% (by volume). For the ideal case of no particle segr egation upon discharge, this plot would be a straight line parallel to the abscissa with an ordinate value of one. Any deviation from this behavior is indicative of segreg ation in the system. The maximu m deviation is obtained for the case of cohesionless (Bo = 0) particles. As the Bond number increases, the discharge segregation decreases. Figure 4(b) presents a similar plot for Bo > 1 and a liquid content of L = 2 %. It is evident from this plot that segregation is negligible at Bo = 1, with larger Bond numbers giving only a minimal change in the segregation profile as compared to the Bo = 1 case. A segregation index may be defined to charac terize the extent of segregation over a full discharge as follows: j(/)1 Segregation index = if jxx N (45) where j is the index indicating the discharge sample and N is the total number of discharge samples. 61 PAGE 62 This segregation index can be norm alized with respect to the segregation index for the cohesionless particle system case as: Segregation index Normalized segregation index = Segregation index for cohesionless case (46) A plot of the normalized segregation index as a function of the Bond number is shown in Fig. 418. There is a decrease in segregation between Bo = 0 and Bo = 1 with the segregation index asymptoting to a minimum value for Bo >1. The error bars represent the 95% confidence interval based on three simulations using differe nt initial conditions. Segregation cannot be totally removed from the system because of the initial nonuniformity of the mixture and sampling error in the data. The primary mechanism responsible for segregat ion in this system is percolation, or sieving, whereby there is a net movement of sma ll particles in the direction of the gravitational acceleration through the interstices be tween larger particles. The ex tent of this mobility of the smaller particles through the matrix of la rge particles depends strongly on cohesion. For the cohesionless case, there exists a flow ing free surface region and discharging core that are finesdepleted. In cont rast, these depletion zones are much more muted in the cohesive case. These trends can be attributed directly to the percolation of fi nes which is much more pronounced in the case of c ohesionless particles. For Bo = 1, the maximum cohesive force on the smaller particles is equal to the weight of the particles. At this value of the Bond number, the smaller particles in the system remain attached to their larger neighbors and percolation ceases. For Bo < 1, the maximum cohesive force on the particle is less than the smaller particles weight and some percolation occurs, re sulting in segregation. For Bo > 1, there is no effect on segregation beyond its critical mini mum since the cohesive force is beyond the critical value that eliminates the gravityinduced percolation. 62 PAGE 63 Effect of Liquid Content on Segregation Liquid bridge force during the enduring contac t is a very weak function of the liquid content. Liquid content, on the other hand, dire ctly affects the range of particle separation distances over which the cohesive force acts. In creasing the liquid content increases the rupture distance and hence, the distance over which cohesion acts in the system. The segregation profile for different liquid contents for a given Bond number ( Bo = 0.5) is shown in Fig. 420. It is important to note that at extremely small liquid contents, some of the assumptions made in the liquid bridge model will become inaccurate. Specifically, the distribution of the liquid on a part icles surface is unlikely to form a thin film that will be available for the formation of a liquid bridge fo r an arbitrary contact point on the particles surface. Since the lower limit at which the distribut ed thin film assumption becomes inaccurate is unknown, an arbitrarily chosen lower value of L = 1% (by volume) was used. The maximum value of L = 23% (by volume) was chosen, which is the upper limit for pendular liquid bridges for a monodisperse system of spherical particles (Yang et al., 2003). These liquid contents are also within the limits for which Mikami et al. (1998) and Soulie et al. (2006) used in their regression expressions. There is a negligible change in the segregation profile as the liquid content, L increases from 1% to 23%. For example, there is only a 4% decrease in the normalized segregation index when the liquid content, L is quadrupled from 4% to 16% Unlike the Bond number, liquid content does not strongly aff ect the segregation profile. Recall that the liquid content governs the rupture distance of the liquid bridge. The flow through the hopper is dense in the sense that particles tend to remain close to each other during discharge, especially in regions far from the ou tlet where percolation takes place. In the hopper system, most of the particles at the onset are in contact with each other a nd as the particles flow 63 PAGE 64 out of the system there is little change in that initial dense configuration. Hence at anytime, except near the hopper outlet, the separation dist ance of a given particle from its closest neighbors remains very small. Thus, even a small amount of liquid is enough to keep the cohesive force sufficiently high to counteract the tendency of the particles to percolate. However, in a dilute granular flow, liquid conten t may play a more important role in particle segregation. Experimental Validation This section presents the results obtained from the experimental investigation and compares them to the simulation results. Specifica lly, three particular characteristics of the flow are compared to the corresponding simulation resu lts: discharge rate, angle of repose of the material left in the hopper and size of clumps formed during the hopper discharge. Discharge Rate Discharge rate is the most important flow char acteristic that helps validate the model. Experiments were performed in both the 90 (flatbottomed hopper) and the 55 hopper. The results for the 90 hopper are summarized in Table 42 and for the 55 hopper in Table 43. A set of three experimental runs were carried out for each system and the 95 % confidence interval from the experiments are listed in the table. In general, there is a very good agreement betw een the discharge rate results obtained from the experiments and the simulations. The simu lation results lie with in the 95 % confidence interval of the experimental result s. One of the important source of error in the experiment is the assumption in the model that the liquid is uniformly distributed. In real systems this is difficult to achieve. At low Bond numbers, as was disc ussed before, the liquid content is not very important in determining the discharge rate. Hence, nonuniformity of liquid content plays a larger role in determining the discharge rate at higher Bond numbers. This nonuniformity of 64 PAGE 65 liquid distribution is the reason for the larg est deviation from the simulation results as is seen in the experiments performed at the highest Bo =3.5. Angle of Repose In the flatbottomed hopper there is always some material left in th e hopper after the flow has stopped. This material left in the hopper ha s a characteristic angle of repose depending on friction coefficient of the materi al and cohesion. Experiments sugge st that the cohesionless glass beads have a characteristic angle of repose sma ller than the cohesionless cast iron shots. This difference can be attributed to the difference in th e particleparticle friction coefficient. The glass particles are simulated with the base case friction coefficient of 0.2, whereas, a higher value of friction coefficient of 0.6 is used to simulate th e cast iron shots. The results from the experiments are compared to the simulation results in Table 44. The simulation results match extremely well with the experimental results. It is observed that the angle of repose increases as the cohesion increases. This me thod can be used as a test to determine the friction coefficients of materials to be used in the DEM simulations. First a simple experiment needs to be run in order to determin e the angle of repose. A parametric study can be done to determine the friction coefficient to be used in the simulati ons that matches the experimental angle of repose. This friction coeffi cient will be the best guess in determining the flow properties of the material. Clump Size The final characteristic of flow investigated in the experimental st udies was the size of clumps formed in the cohesive cases. For the case of cohesionless particles, no clump formation was seen in both the simulations and experiments. Figure 421 s hows two typical clumps formed in three different cases of Bond number corr esponding to the four experimental systems described in Table 31. 65 PAGE 66 It is observed that the size of clum ps increa ses as cohesion increases. On an average, it is observed that for Bo =0.3, maximum clump sizes of 34 particles forms. For Bo =1, maximum clump sizes of the order of 1516 partic les are seen, where as increasing the Bo=3.5, increases the maximum clump sizes to the order of 35 particle s. It should be noted that smaller clumps are also present in the discharge for Bo=3.5, but the maximum clump size in the distribution is of the order of 35 particles. The experimental results obtai ned from a high speed camera at 400 fps are shown in the figure 422. The size of clumps can be calculated by using pixel anal ysis described in chapter 3. The results match well with the simulation results for all the four cases. For Bo =0.3 shown in figure 422A, for the cast iron shot s with silicone oil, a maximu m clump size of 34 particles are observed as in the simulations. Both the Bo =1 cases, i.e., glass beads with silicone oil or cast iron shots with water, yield a maximum averag e clump size of 14 particles. For the highly cohesive case of Bo=3.5 corresponding to the gl ass beads with water, a maximum average clump size of 35 particles is observe d as in the simulations. 66 PAGE 67 A Figure 41. Discharge profile for different friction coefficients. A) = 0.05 and w = 0.05 and B) = 0.84 and w = 0.84. The value in parentheses is the mass discharge rate. 67 PAGE 68 B 68 PAGE 69 A Figure 42. Instantaneous mass discharge rate plotted as a function of instantaneous dimensionless fill height for different friction coefficients. A) = w = 0.05.B) = w = 0.1 C) = w = 0.15 D) = w = 0.2 E) = w = 0.5 F) = w = 0.84 69 PAGE 70 B Figure 42. Continued 70 PAGE 71 C Figure 42. Continued 71 PAGE 72 D Figure 42. Continued 72 PAGE 73 E Figure 42. Continued 73 PAGE 74 F 74 PAGE 75 Figure 43. Mass discharged from the hopper plotted against time fo r a range of particleparticle and particlewall friction values. The value in parentheses is the mass discharge rate. 75 PAGE 76 Figure 44. Mass discharged from the hopper plotted against time for two values of the coefficient of restitution. 76 PAGE 77 Figure 45. Mass discharged from the hopper plotted against time for hoppers with two different widths, W 77 PAGE 78 Figure 46. DEM mass discharge ra te plotted against hopper angl e (from the vertical). The correlations by Brown and Richar ds [Eq. 14] and Rose and Tanaka [Eq. 15] are also plotted. 78 PAGE 79 Figure 47. Mass discharg e rate per unit length raised to the 2/3 power plotted as a function of the outlet width, W0. The solid line is the trend pr edicted by the Beverloo correlation. 79 PAGE 80 Figure 48. Mass discharge from the hopper plot ted against time for tw o different particle diameters, d 80 PAGE 81 Figure 49. DEM mass discharge ra te plotted against the fines mass fraction for a bidisperse system. The correlation of Humby et al [Eq.17] is also shown. 81 PAGE 82 Figure 410. Dimensionless mass di scharge rate as a function of Bond number for various values of dimensionless outlet width and L = 2% 82 PAGE 83 Figure 411. Mass discharge rate per unit length raised to the 2/3 power plotted as a function of the outlet width, W0, for various values of Bond number and L = 2 % 83 PAGE 84 Figure 412. Beverloo constant, k as a function of Bond number for L = 2% 84 PAGE 85 A Figure 413.Snapshot of disc harge. A) Cohesionless, Bo = 0, Wo/d =10.6, L = 2% B) Cohesive, Bo = 2.5, Wo/d =10.6, L = 2 % 85 PAGE 86 B 86 PAGE 87 Figure 414. A comparison of DE M simulation results with predictions from the modified Beverloo correlation for cohesion 87 PAGE 88 A Figure 415. Mass discharged from the hopper as a function for time for different Bond Number. The values in parentheses show the ma ss discharge rate per unit length. A) Bo = 0.3 B) Bo = 2 88 PAGE 89 B 89 PAGE 90 Figure 416. Mass discharge rate plotted agains t hopper angle (from the ve rtical) for different Bond numbers with Wo/d = 10.6 and L = 2 %. 90 PAGE 91 A Figure 417. Discharge segrega tion for varying Bond numbers. L = 2%, f x = 5%, D = 4.3, and = 90 A) for Bo < 1 B) for Bo > 1. 91 PAGE 92 B 92 PAGE 93 Figure 418. The normalized segregation index plotted as a function of Bond number. 93 PAGE 94 A B Figure 419. Local fines fraction fields. A) Bo = 0 (cohesionless), and B) Bo = 1. Red indicates a large fines fraction as compared to a well mixed state, and blue indicates a small fines fraction. Gray regions indicate well mixed regions. 94 PAGE 95 Figure 420. Discharge segrega tion for varying liquid content ( L ), with Bo = 0.5, f x = 5%, D = 4.3, and = 90. 95 PAGE 96 A B C Figure 421. Illustration of Clump sizes from simulation. A) Bo=0.3 B) Bo =1 C) Bo =3.5 96 PAGE 97 A B Figure 422. Illustration of Clump sizes from experiments. A) Silicone oil with Cast iron, Bo =0.3 B) Water with Cast iron, Bo = 1 C) Silicone Oil with Glass, Bo =1 D) Water with Glass, Bo =3.5 97 PAGE 98 D Table 41. List of various frictional cases Case w 1. Baseline 0.2 0.2 2. Large particlewall friction 0.2 0.84 3. Large particleparticle friction 0.84 0.2 4. Large particle and wall friction 0.84 0.84 Table 42 Comparison of experimental discharg e rate to that of simulations in the 90 hopper System Experimental Average Discharge Rate 95 % Confidence Interval Simulation Results Glass Beads Bo = 0 139.4 g/cm/s 16.7 140.4 g/cm/s Cast iron Shots Bo=0 394.6 g/cm/s 89.8 459.6 g/cm/s Cast iron Shots with Silicone Oil Bo = 0.3 381.78 g/cm/s 53.1 436.44 g/cm/s Glass Beads with Silicone Oil Bo = 1 115.5 g/cm/s 7.74 118 g/cm/s Cast iron Shots with Water Bo = 1 350.43 g/cm/s 61.70 408 g/cm/s Glass Beads with Water Bo = 3.5 71.9 g/cm/s 19.8 24.1 g/cm/s 98 PAGE 99 99 Table 43 Comparison of experimental discharg e rate to that of simulations in the 55 hopper System Experimental Average Discharge Rate 95 % Confidence Interval Simulation Results Glass Beads Bo = 0 140.7 g/cm/s 14.4 138.2 g/cm/s Cast iron Shots Bo=0 416.93 g/cm/s 75.3 454.5 g/cm/s Cast iron Shots with Silicone Oil Bo = 0.3 393.72 g/cm/s 51.5 442.4 g/cm/s Glass Beads with Silicone Oil Bo = 1 121.86 g/cm/s 6.24 125.8 g/cm/s Cast iron Shots with Water Bo = 1 347.88 g/cm/s 72.3 408.6 g/cm/s Glass Beads with Water Bo = 3.5 75.52 g/cm/s 21.2 52.64 g/cm/s Table 44. Comparison of experi mental results for angle of re pose to that of simulations System Experimental Angle of Repose Simulation Results Glass Beads Bo = 0 17 19.65 Cast iron Shots Bo=0 26.5 26.53 Cast iron Shots with Silicone Oil Bo = 0.3 28 29.7 Glass Beads with Silicone Oil Bo = 1 31.5 32.1 Cast iron Shots with Water Bo = 1 37.8 36.70 Glass Beads with Water Bo = 3.5 39.80 47.8 PAGE 100 CHAP TER 5 CONCLUSION AND RECOMMENDATIONS This chapter discusses the key conclusions draw n from the results presented in Chapter 4. Then, it discusses the avenues for furt her research in the next section. Conclusions The study of cohesionless discharge rate uses the discrete element method to model the discharge of bulk granular materi als from a quasithreedimensiona l, wedgeshaped hopper. The effects of various hopper and particle properties on discharge rate were studied and the results were compared to experimental correlations report ed in the literature. In general, there is good agreement between the DEM disc harge rate predictions and th e values from experimental correlations for all of the inve stigated parameters. These findi ngs provide strong evidence that DEM simulations can be used to predict the flow behavior of granular materials through hoppers. The coefficient of restitution was shown to have negligible influence on the discharge rate. In contrast, friction strongly affects the discharge rate. In particular, particleparticle friction is more important in determining the discharge rate than particlewall friction. For friction coefficients greater than 0.2, the discharge rate is independent of fill height for the systems investigated here. The DEM simulations also sh owed that the discharge rate varies with exit width (corrected for the dead space, Z = kd ) raised to the 3/2 power, identical to the predictions of the modified Beverloo equation for rectangul ar hoppers. The influence of hopper angle also matched well quantitatively with the correlations of Brown and Richards (1965) and Rose and Tanaka (1959), with the Rose and Tanaka ( 1959) correlation matching the DEM results more closely at smaller angles and the Brown and Ri chards (1965) correlation working better above a critical angle. Lastly, the discharge rates obtained from DEM simulations for bidisperse systems 100 PAGE 101 agree well with th e experimental correlations of Humby et al. (1998), although the correlation gives values that are consistently larger than the DEM values. In terms of practical applica tions, the results indicate that hopper discharge rate can be increased by decreasing particle friction, as we ll as increasing outlet wi dth. Hopper angle can also be used to increase the discharge rate from a given hopper geometry (for a fixed outlet width), provided the chosen hopper angle is below the critical angle. In a ddition, it is shown that in order to maintain a uniform discharge rate, th e fill height should not fa ll below one half of the outlet width. Finally, the simula tions confirm that particle si ze hopper width, coefficient of restitution, and particlewall friction (assuming a fl at bottomed bin) have a negligible effect on hopper discharge rate. The study of cohesive discharge rate employs a pendular liquid bridge force to model cohesion between particles. The effects of vary ing Bond number, outlet width, liquid content and hopper angle on the discharge rate ar e studied. The results show that the basic discharge rate trend given in the Beverloo correlation, where the mass discharge rate varies linearly with the (corrected outlet width)3/2 is also valid in the case of cohesive particles. Using the method used by Beverloo et al. to develop their correlation, a similar correlation was developed for the case of cohesive discharge with a modified Beverloo constant which is a function of Bond number. The prediction from this modified correlation provides a reasonable match to the simulation data, with more accurate estimates at low to moderate Bond numbers. This correlation can likely be improved with a better estimate of the flowing bulk density as a function of cohesion. The influence of liquid content, which st rongly affects the distance over which liquid bridge forces act, was found to be dependent on Bond number. At low Bond number, identical 101 PAGE 102 discharg e rates were obtained despite varying liquid content over a wide range. However, for high Bond number, there was a slight decrease in the discharge rate with increasing liquid content. Finally, for cohesive material, the in fluence of hopper angle on discharge rate followed a trend very similar to the c ohesionless case, although the critical flow angle decreased as cohesion increased. The study of cohesive segregation uses the discrete element method to model the discharge of cohesive bidispersed gra nular materials from a quasithreedimensional, wedgeshaped hopper. A pendular liquid bridge fo rce is used to model cohesion between particles. The effects of varying Bond number and liqui d content on the extent of disc harge segregation are studied. For Bond numbers equal to or greater than one, with the Bond number based on the smaller particles, the maximum cohesive force on a smaller particle is greater than or equal to the particles weight. As a result gravityinduced percolation of the small particles through the larger particles is sign ificantly reduced and discharge segr egation from the hopper is minimized. A Bond number greater than one has little additional effect on segregation since the cohesive force acting on a small particle is already suffi ciently large to overcome the small particles weight. However, Bond numbers greater than on e do continue to have a significant influence on the discharge rate from the hopper, with the fl ow decreasing within in creasing Bond number. Hence, a hopper flow with a Bond number of one appears to be the optimal case which minimizes segregation while at the same time producing a discharge rate that is larger and steadier than those for Bond numbers greater than one. Finally, because of the dense flow pattern existing in a hopper system, it was found that liquid content played a negligible role in determining the extent of discharge segregation. 102 PAGE 103 Finally, experim ents are performed in hoppers of identical dimensions as in the simulation to verify the model. Three specific characteristic of the flow: discharge rate, angle of repose of the material and size of clumps formed during th e discharge are investigat ed experimentally and compared to the simulations. The simulation results agree well with the experimental results. The angle of repose calculated from the experiments can be used in a parametric study, to get a physically more realistic frictiona l coefficient to be used with the DEM simulations for real systems. In general, the discha rge rate decreases, angle of repose increases and size of clumps increases as cohesion increases. Recommendations The present work has reached some important conclusions regarding hopper discharge rate and segregation and the parameters that affect it but there are several areas open to further research. These include extending the computational model for nons pherical particle shapes and using a different cohesion model. Each of this is discussed. The current DEM model is valid only for spherica l particles. In real situations present in the industry many of the systems are nonspherical in nature. By incorpor ating particle shape effects into the model, its applicability can be significantly broadened to include a variety of industrially relevant systems. Some commonly used nonspheri cal shapes in 2D include continuously connected disks (Popatov and Campbell, 1998), polygons (Matuttis et al ., 2000), and ellipses (Dziugys and Peters, 2001) and, in 3D, rigid clusters of spheres in various shapes (Walton and Braun, 1993), superquartics (Williams and Pentland, 1992), dilated disks (Hopkins and Tuhkuri, 1999) and ellipsoids (Lin and Ng, 1997). As contact detection typically becomes computationally intensive for nonsp herical shapes, this should be a key criterion when selecting a nonspherical shape. As such, rigid clusters of spheres may be the most attractive, since contact detection is still based on spherical shapes. 103 PAGE 104 104 The pendular liquid bridge model used fo r cohesion in this work is the most computationally intensive method which offers th e advantage that all the required parameters (surface tension, liquid content a nd contact angle) are experiment ally measurable quantitities. However, this model has limitations of its own. It is valid only for a range of liquid contents. At extremely low liquid content, the model breaks down because the assumption of formation of distributed thin film over the particle becomes inaccurate. At high liquid content, the system moves into the funicular regime where the liqui d bridge interactions are no longer binary. Several implementations of cohesive model that are relatively simple such as those using a cohesion constant similar to a spring stiffness [Asmar et al., 2002 and Luding et al., 2003] or a cohesive force proportional to the contact length (Matuttis an d Schinner,2001) exist in the literature. Others have also used a square well (Weber et al ., 2004) or continuous interaction potential with an attractive well (Baxter et al ., 2000). A comparison of implementation of one of these simpler models with the more elaborate pe ndular liquid bridge m odel can give a deeper insight into the cohe sive flow phenomenon. It is also possible to do a more accurate and comprehensive study of clump sizes formed during discharge both using experime nts and the simulations. The basic trends seen in the results shown in the largest clump sizes formed during discharge suggests a corr elation with the Bond number. PAGE 105 APPENDIX A SAMPLE I NPUT FILE 6790 np: total number of particles/ 1 ifill: flag: 1=random, 2=hard coded, 3=dense grid, 4=file input, 5=RHCP, 7=FoverC, 1 randit: selects various random configurations, if ifill=1/ 1 nreps: number of replicates (with different random ICs), 2 ibc: flag: 2=Hopper (plane flow), 3=Hopper (cylindrical), 4=Hopper (pharma) / 0 nrecyc: number of particles to recycle (enter 2 for 2*np), 0=no recycle/ 1 irot: flag: 1=initial random rotation on, 0= off/ 3 isearch: flag: 1=brute force, 2=grid, 3=radius(NN)/ 3 idim: flag: dimensionality: 2D or 3D/ .56 z_width: width of zdirection periodic boundaries for 3D hopper (ibc=2)/ 22.5 uho: upper hopper origin upper hopper variables used only for ifill=9/ 21.5 hht: upper hopper height/ 10 hdia: upper hopper diameter (width)/ 1.3 hod: upper hopper orifice diameter/ 15 hha: upper hopper halfangle, from vertical, in degrees/ 30 hht: hopper height/ 12.5 hdia: hopper diameter (width)/ 2.5 hod: hopper orifice diameter/ 89.99 hha: hopper halfangle, from vertical, in degrees/ 1 idiam: flag: particle diameters: 1=monodisperse, 2=bidisperse/ 2 idd: flag: particle diameter dist'n: 1=uniform, 2=Gaussian, 3=lognormal/ 0.235 dia: diameter of monosized particles or larger particles, if binary system/ .01 sd: std dev. of monosized particles or larger particles, if binary system, and idd>1/ 7.85 dens: density of monosized particles or larger particles, if binary system/ 6790 nLG: number of large particles/ 0.0521 dia_sm: diameter of smaller particles, if binary system/ 0.0048 sd_sm: std dev. of smaller particles, if binary system, and idd>1/ 7.85 dens_sm: density of smaller particles, if binary system/ 80 num_sam: number of "samples" to take at initial fill/ 0.5 sam_rad: radius of sample/ 0.00001 dt: initial time step/ 0 itadj: flag: 1=adjust dt wrt to overlap, 0=constant dt/ 5000000 tfinal: number of time steps/ 2500 twrite: interval for writing to file/ coh.mono name: file name prefix for data output/ 0.2 mu_ll: coefficient of fr iction: largelarge contacts/ 0.2 mu_ss: coefficient of fr iction: smallsmall contacts/ 0.2 mu_sl: coefficient of fr iction: smalllarge contacts/ 0.01 mu_r: rolling coefficient of friction/ 0.2 mu_lw: pw coefficient of friction: largewall contacts/ 0.2 mu_sw: pw coefficient of friction: smallwall contacts/ 0.01 mu_rw: pw rolling coefficient of friction/ 250000 k1_ll: loading spring constant: largelarge contacts/ 105 PAGE 106 283000 k2_ll: unloading spring constant: largelarge contacts/ 250000 k1_ss: loading spring constant: smallsmall contacts/ 283000 k2_ss: unloading spring constant: smallsmall contacts/ 250000 k1_ls: loading spring constant: largesmall contacts/ 283000 k2_ls: unloading spring constant: largesmall contacts/ 0.2857 kt_r: kt to k1 ratio for pp contacts/ 250000 k1_lw: wall loading spring constant: walllarge contacts/ 308600 k2_lw: wall unloading spring constant: walllarge contacts/ 250000 k1_sw: wall loading spring constant: wallsmall contacts/ 308600 k2_sw: wall unloading spring constant: wallsmall contacts/ 0.2857 ktw_r: ktw to k1w ratio for pw contacts/ 1 ivdamp: flag: 1=turns on viscous damping during settling phase, 0=off/ 1 igrav: flag: 1=gravity on (ydir), 0=gravity off/ 981 gravacc: gravitational acceleration (cm per sec per sec)/ 0 ivdist: flag: 1=print velocity distributions, 0=off/ 0 ifile: flag: 1=print to file(s), 0=off/ 1 imovie: flag: 1=write movie/ 1 ihe: flag: 1=output hopper energies to screen/ 1 isup: flag: 1=suppress all other output, 0=print other/ 0 ipos: flag: 1=output positionvelocity, 0=off/ 0 icoll: flag: 1=output collision details on, 0=off/ 72 gamma_w: Surface Tension of liquid/ 0.0001 viscosity_w: Viscosity of water/ 0.1 theta_wg: Contact angle waterglass/ 1.256 theta_ww: Contact angle waterwall/ 0.00082 Volume1: Volume of water with large bead/ 0.00082 Volume2: Volume of water with smaller bead/ 0.00224 asperity: one hunredth of particle diameter/ 0.0 Volume_Wall: Volume of water supplied by the wall for Liquid Bridges/ 106 PAGE 107 APPENDIX B HOPPE R DISCHARGE CODE /* 3D HOPPER FLOW of Bidisperse spheres in a Wedge Hopper with cohesion using Liquid Bridges.*/ #include PAGE 108 // INDEX VARIABLES int b,d,c,h,i,j,k,p,m,u,w; int bb, ir; double aa, cc, dd, kk, mm, nn, pp, qq, rr, ss, tt, uu, vv, ww; int rand_factor, randit, nreps; int curr_randit=1; int indy; // SYSTEM INPUT VARIABLES int np, nLG; double npinv, dtinv; int ifill, ibc, irot, isearch, idim, idiam, idd, itadj, ivda mp, igrav, ipos, ivdist, ifile, imovie, ihe, isup, icoll, iopen; int irf; double dia, dia_sm, sd, sd_sm, dens, dens_sm, dia_ln, dia_sm_ln, sd_ln, sd_sm_ln; double xmax, ymax; // PARTICLE PROPERTIES/LOCATIONS/VELOCITIES double *x, *y, *z, *cx, *cy, *cz, *wx, *wy, *wz, *eta, *xi, *lambda, *mass, *diam, *moment, *w_magn; //double *tempcx, *tempcy, *tempcz, *oldtempcx, *oldtempcy, *oldtempcz; double tempz; double wi_magn, wj_magn; // magnitude of rotational velocities double wx_i, wy_i, wz_i; // rotational velocities in the inertial frame double wx_b, wy_b, wz_b; // rotational velocities in the bodyfixed frame double totalmass = 0.0; //total mass of hopper load double mindiam=10.0; //record min diam double maxdiam=0.0; double maxmass=0.0; double minmass=100.0; double avediam, avediamSM, avediamLG; // CONTACT SEARCH VARIABLES int *pw_nn_list, *skip_part, *skip_ search, *pp_collisions, *pw_collisions; double** separtn; double xz_sq, xz_off_sq, xz_varoff_sq, par; double pw_max_dist, pw_maxvel, test_vel; int pw_neighbors; int pwn_updates=0; double search_count[5];// needs to be large in itially so neighbors are determined at t=1 !! int n_updates=0; int nw_updates = 0; int issmall=0; int islarge=0; int issmall2=0; int islarge2=0; int startA, endA, startB, endB; // CONTACT COUNTER VARIABLES int wall_collisions[9]; // pw collision counter (spec'd by wall) int n_col = 0; // position in separtn[][] "table" int part_collisions = 0; int part_collisions_set = 0; int total_wall_collisions = 0; 108 PAGE 109 int total_wall_collisions_set = 0; int pp_collisions_segn_stats = 0; int pw_collisions_segn_stats = 0; // OVERLAP STATISTICAL VARIABLES double overlap=0; double min_overlap=1; double max_overlap=0; double sum_overlap=0; double min_pwoverlap=1; double max_pwoverlap=0; double sum_pwoverlap=0; double min_overlap_set, max_overlap_set, sum_overlap_set, min_pwoverlap_set, max_pwoverlap_set, sum_pwoverlap_set; // HARDPARTICLE CONTACT MODEL VARIABLES double g[4]; // relative velocity at contact point double n[4]; // normal unit vector for contacts double s[4]; // tangential unit vector for a given collision. double n_magn, s_magn, g_dot_n, g_dot_s, nm, sm; double fcn_overlap; // overlap given by trendline fit of guidelines from Hopkins & Louge 1991 // SOFTPARTICLE CONTACT MODEL VARIABLES double Fnij[4], Ftij[4], torqueij_i[4], torqueij_b[4]; // normal force, tangential fo rce, inertial frame torque, and bodyfixed frame torque double Fnji[4], Ftji[4], torqueji_i[4], torqueji_b[4]; double** Fnsum; // sum of forces acting on a given particle double** torquesum; // sum of torques acting on a given particle int** neighbors; // number (in 1st column) and list of neighbors of a given particle double Fn_magn, Ft_magn; double delta_o; double delta; double tan_dispx, tan_dispy, tan_dispz; double psi_total, psi_x, psi_y, psi_z, psi_magn; int num_integrations=0; int pp_total_integrations=0; int pw_total_integrations=0; double ovlp_for_se=0.0; int countlargedelta=0; int overlap_warning=0; int overlap_warning_set=0; // CONTACT COEFFICIENT VARIABLES double mu, muw; double mu_ss; // Friction for smallsmall contacts double mu_ll; // Friction for largelarge contacts double mu_sl; // Friction for smalllarge contacts double mu_r; // Rolling resistance for pp contacts double mu_sw; // Friction for smallwall contacts double mu_lw; // Friction for largewall contacts double mu_rw; // Rolling resistance for pw contacts double kt, ktw; 109 PAGE 110 double k1, k2, k1w, k2w, k1_ll, k2_ll, k1_ss, k2_ss, k1_sl, k2_sl, k1_lw, k2_lw, k1_sw, k2_sw; // Spring constants; double k_ratio, kt_f, ktw_f; // TIME STEP VARIABLES double dt; // Time step size int t=0; // index: # of time steps ==> time initial time = t dt (for constant dt) int tfinal; // number of time steps int twrite; // interval to write particle state info double total_time = 0.0; // Total Real time of simulation (seconds) int topen=9999800; // step at which to open hopper gate (after a quasistatic state has been reached). Set to large number initially double hrs, mins, secs; double clocktime; // HOPPER GEOMETRICAL VARIABLES double uho, uhht, uhdia, uhod, uhha, uhhar, hht, hdia, hod, hha, hhar, sinhhar, coshhar, tanhhar, tanuhhar; double sinb, cosb, sinbp1, cosbp1; double h1, h2, alph, gam, Hdist; //double y_shift = 0.0; // shift for pharma hopper to account for slide gate within discharge chute (11cm above opening) // HOPPER PERIODIC BOUNDARY VARIABLES double z_width; // total width of zdir Periodic Boundaries (from z=0.5*z_width to z=0.5*z_width) int corz; // LEESEDWARDS PBC VARIABLES double per_disp, total_disp; // displacem ent values for "moving" periodic images int cory, corx, coryy, corx2; // indices to determine if particle is crossing top vs. bottom, L vs. R double temp2; // temporary variables for storage of x[i] and cx[i] for use with LEPBC int greg = 0; // SYSTEM ENERGY VARIABLES double ve, ke, re, te, pe, se; //kinetic, rotational, total, potential, and spring energies // STRESS TENSOR VARIABLES //double Tcave[4][4], Tkave[4][4]; //double fluc[4][4]; //double cxave, cyave, avecz, czave; //double aveT[4][4], aveTc[4][4], aveTk[4][4]; //double flucx = 0, flucy = 0, flucxy = 0; //double instJn[4][4];; //double perchangeTk[4][4]; //double perchangeTc[4][4]; //double maxperchange=0.0; //double instT[4][4], instTc[4][4], instTk[4][4]; //double instaveT[4][4], instav eTc[4][4], instaveTk[4][4]; //double instfluc[4][4]; //double instflucx = 0, instflucy = 0, instflucxy = 0; //double previnstaveT[4][4], previnst aveTc[4][4], previnstaveTk[4][4]; //double Jn[4][4], Jntot[4][4]; 110 PAGE 111 //double temperature; // Granular Temperature //double insttemperature=0.0; // PP OVERLAP DISTRIBUTION VARIABLES int num_bin_o=200; double deltabin_o; double bin_o[301][3]; // PW OVERLAP DISTRIBUTION VARIABLES int num_bin_pwo=200; double deltabin_pwo; double bin_pwo[301][3]; // PARTICLE SIZE DISTRIBUTION VARIABLES int num_bin_psd=50; double deltabin_psd; double bin_psd[61][3]; int sum_of_psd = 0; // AVERAGE SIFTING VELOCITY VARIABLES int num_bin_x_sift=0; int num_bin_y_sift=0; double deltabin_sift; double bin_sift[160][140][13]; double c_x, f_x, c_y, f_y; // HOPPER DISCHARGE POSITION DISTRIBUTION VARIABLES int num_bin_hdp=15; double deltabin_hdp; double bin_hdp[60][7]; int sum_of_hdp = 0; // HOPPER DISCHARGE AVERAGING VARIABLES "BY MASS" DISTRIBUTIONS int num_disch_in_massbin = 0; int distnlocation = 1; double flowdistn[201][15][102]; double avedischdiam, avedischvel, avedischmass, avedischmassS, avedischmassL, avedischdens; double cummassdisch=0.0; double expecteddiam=0.0; double expectedmass=0.0; int mass_update_time, old_mass_update_time; // HOPPER DISCHARGE AVERAGING VARIABLES "BY TIME" DISTRIBUTIONS int num_disch_in_timebin = 0; int tdistnlocation = 1; double tflowdistn[201][15][102]; double tavedischdiam, tavedischvel, tavedischmass, tavedischmassS, tavedischmassL, tavedischdens; double tcummassdisch=0.0; int avedischtime = 0; // APPROACH VELOCITY DISTRIBUTION VARIABLES int num_bin_av1=50; double deltabin_av1; double bin_av1[151][3]; int num_bin_av2=50; double deltabin_av2; 111 PAGE 112 double bin_av2[151][3]; // VELOCITY DISTRIBUTION VARIABLES int num_binf=22; int num_bina=20; double deltabinf, deltabina; double binfx[125][3], binfy[125][3]; double binax[125][3], binay[125][3]; // RHCP FILL VARIABLES double xinc, yinc, old_yinc, spacing; // OUTPUT FILES VARIABLES string movname; string exec_name; string filename1, filename2; string filename3; string filename4, filename5; string filename6; string filename7; string filename8; string filename9; string filename10; string filename12, filename13; string filename16, filename17, filename18; // track particle velocity distributions string filename21, filename22, filename23, filename24, filename25, filename26, filename27; string filename29, filename30, filename31; string filename32, filename33,filename34,filename35,filename36; string zero="0"; string doublezero="00"; // PARTICLE VELOCITY TRACKING VARIABLES double maxrot=0.0; int fastest_rot_p=1; double maxvel, maxxvel, maxyvel, maxzvel=10; double oldmaxvel=10; // Maximum particle velocity used for calculation of adjusted dt double oldmaxvel_L=10; double oldmaxvel_S=10; int fastest_p, fastest_p_L, fastest_p_S, fastest_x_p, fastest_y_p, fastest_z_p; // PARTICLE RESIDENCE TIME DIST'N VARIABLES double t_star, time_k; double mindischtime=10.0; double maxdischtime=0.0; // vars for hopper RTD colormap double fixed_time; // AVERAGE VELOCITY COLORMAP VARIABLES double v_star; double *ave_velocity; int frame2=0; int frame3=0; int num_velocities=1; // don't start at 0, else NaN error for first frame //Anshus //LIQUID BRIDGE COLORMAP VARIABLE 112 PAGE 113 int frame4=0; //Anshue // "COUNT" VARIABLES int writecount6 = 0; int writecount8=15000000; int count3 = 0; int count8 = 0; // VISCOUS DAMPING VARIABLES slow particles durin g freefall that occurs during filling to reduce overall runtime double vdc = 0.999; // parameters for viscous damping during particle settling phase double initial_pe; // HOPPER FLOW RECYCLING VARIABLES recycle discharged particle s to reach a steady state flow rate double nrecyc; // input variable that define s number of particles to recycle listed as a multiple // of 'np'... e.g. nrecyc=2 im plies that 2*np particles will be recycled int num_recycled = 0; // a counter wh ich indicates the number of particles that have currently been recycled int is_recycled = 0; // flag that a particle has been recycled; update of neighbor lists is required int done = 0; int contact = 0; double recyc_dist; double** discharge; // PW OVERLAP VARIABLES int wall_num; double** pwdeltao; // for pw collisionswalls are numbered from 0 to 8 in a CW manner starting w right bin wall: double** pwtan_dispx; double** pwtan_dispy; double** pwtan_dispz; /* Wall Wall Number ========================================= Right Bin Wall 0 Right Hopper Wall 1 Pt btn RHW and ROW 2 Right Orifice Wall 3 Slide Gate 4 Left Orifice Wall 5 Pt btn LOW and LHW 6 Left Hopper Wall 7 Left Bin Wall 8 */ // "OTHER" VARIABLES inline double sq(double s) {return s*s;} inline double cube(double s) {return s*s*s;} double sarea, varea, sfrac, vfrac, svol, vvol, totalvol, totalarea; 113 PAGE 114 double Pi = 2*acos(0.0); double t1, t2, delta_time; double eta_dot, xi_dot, lambda_dot; double floorht=150; int num_discharge=0; int prev_disch_time=0; double x_fines; int maxneighbors=0; double ave_neighbors=0; double maxdelta=0; int loopstop, loopstop2; double grav; // Grav accel to use in calculations (as opposed to input "gravacc") double gravacc; // input gravitational acceleration double ind; double lt = 0.01; // accounts for 1/2 thickness of line displayed in visualization double y_loc, diamA, diamB; // vars for layered fill int frame; // for pvis colormap double xiprofile[22][4]; double cur_sep; int *discharget; //records time which each particle exits hopper const int large_num = 99999999; /*FUNCTION DEFINITIONS */ //Anshus void ParticleWallLiquidBridgeModel(int i,int wall_num,double _hdia,double _hod, double _hhar); //Anshue //Anshus void RemoveLiquidBridges(int i,LiquidBridge*l); //Anshue void RemoveCollision(int i, collision *a); double HRand(double n) // returns random numbers on [1,1] { return double( 2.0 (1.0 rand() / RAND_MAX) 1.0 ); } double GRand(double mu, double sd) // returns a random number from a Gaussian dist'n with average, mu, and standard deviation, sd { double ff,hh,gg; 114 PAGE 115 do { ff = HRand(0.1); hh = HRand(0.1); gg = 1.0 ff ff + 1.0 hh hh; } while ( gg >= 1.0); gg = sqrt( 2.0 log(gg) / gg); return 1.0 sd ff gg + mu; } string itos (int integer) { stringstream ss; ss << integer; return ss.str(); } double min(double a, double b) { if (ab) return a; else return b; } int sign(double a) { if (a<0) return 1; else return 1; } collision *InsertCollision(int ii, int jj) { int k; collision *a;//, *b; if(ii>jj){ k = ii; ii = jj; jj = k; } a = collision_map[ii]; if(!a){ 115 PAGE 116 a = (collision*)malloc(sizeof(collision)); collision_map[ii] = a; a>parent = 0; a>next = 0; a>j = jj; a>delta_o = 0.0; a>tandx = a>tandy = a>tandz = 0.0; } else{ while(a>j != jj && a>next) a=a>next; if(a>j!=jj){ a>next = (collision*)malloc(sizeof(collision)); a>next>parent = a; a = a>next; a>next = 0; a>j = jj; a>delta_o = 0.0; a>tandx = a>tandy = a>tandz = 0.0; } } a>time = t; //cout<<"TIme STEP COLL "< PAGE 117 l>nex>par = l; l = l>nex; l>nex = 0; l>j = jj; } } //cout<<"TIme STEP in insert Liquid Bridge function is "< PAGE 118 collision *a; // printf("\n"); for(ii=1;ii<=np;ii++) for(a=collision_map[ii];a;a=a>next) { if(a>time==t) printf("%4d %4d %4d %e %e %e %e %e\n",ii,a>j,t,a>cur_se p,a>delta_o,a>tandx,a>tandy,a>tandz); else RemoveCollision(ii,a); } } void ResetWallData() { x[0]=y[0]=z[0] = 0.0; wx[0] = wy[0] = wz[0] = 0.0; // define the wall 'properties'; store in the zero index cx[0] = cy[0] = cz[0] = 0.0; diam[0] = 100000.0; mass[0] = diam[0]*diam[0]; moment[0] = 900000000.0; } void CalcTib(int i) { /* Tib[i][0][0] = cos(xi[i])*cos(lambda[i]); Tib[i][0][1] = cos(eta[i])*sin(lambda[i]) + sin(eta[i])*sin(xi[i])*cos(lambda[i]) ; Tib[i][0][2] = sin(eta[i])*sin(lambda[i]) cos(eta[i])*sin(xi[i])*cos(lambda[i]) ; Tib[i][1][0] = 1.0*cos(xi[i])*sin(lambda[i]) ; Tib[i][1][1] = cos(eta[i])*cos(lambda[i]) sin(eta[i])*sin(xi[i])*sin(lambda[i]); Tib[i][1][2] = sin(eta[i])*cos(lambda[i]) + cos(eta[i])*sin(xi[i])*sin(lambda[i]); Tib[i][2][0] = sin(xi[i]); Tib[i][2][1] = 1.0*sin(eta[i])*cos(xi[i]); Tib[i][2][2] = cos(eta[i])*cos(xi[i]); */ //cout < PAGE 119 } } void RandomFill(double y_origin, double _hht, double _hdia, double _hod, double _hha) { double _hhar, _tanhhar, _sinhhar; _hhar = _hha Pi / 180; _sinhhar = sin(_hhar); _tanhhar = tan(_hhar); cout.setf(ios::showpoint); cout <<"\n starting Random Fill( y_orign, _hht, _hdia, _hod, _hha = "< PAGE 120 if (idim == 3 && ibc==2) { z[i] = HRand(rand_factor)*1.0*z_width; } //2.0*z_width 1.0*z_width;} if (idim == 3 && ibc==3) { z[i] = HRand(rand_factor)*0.55*_hdia; xz_sq = sqrt( sq(x[i]) + sq(z[i]) ); } // cout << "x,y,z = "< PAGE 121 _tanhhar = tan(_hhar); alph = 0.5*(_hdia _hod); cout.setf(ios::showpoint); cout <<"\n starting Random Fill( y_orign, _hht, _hdia, _hod, _hha = "< 0.5*_hdia0.5*diam[1]) ); */ x[1] = 0.0; y[1] = 1.0; z[1] = 0.0; cout << endl<<1; //place particles 2 through 'np' for (i=2; i<=np; i++) { if (np>800 && i%100==0) cout << endl<1.5*_hht  y[i]PAGE 122 } }//close for j }//close for i } void FillLayer(int start, int end, double diameter) { int ii, jj; for ( ii=start; ii<=end; ii++) { contact = 0; do { count3++; if (ibc==2) { do x[ii] = HRand(0.1) 0.55 hdia; while ( fabs(x[ii]) > 0.5 hdia 0.5 diam[ii]  fabs(x[ii]) > tanhhar*y_loc + 0.5*hod 0.5*diam[ii]/coshhar ) ; y[ii] = y_loc; do z[ii] = HRand(0.1) 0.55 z_width; while ( fabs(z[ii]) > 0.5 z_width ); } if (ibc==3) { do { x[ii] = HRand(0.1) 0.55 hdia; z[ii] = HRand(0.1) 0.55 hdia; xz_sq = sqrt( sq(x[ii]) + sq(z[ii]) ); } while ( xz_sq > 0.5 hdia 0.5*diam[ii]  xz_sq > tanhhar*y_loc + 0.5*hod 0.5*diam[ii]/coshhar ); y[ii] = y_loc; } if (ii>1) { for (jj=1; jj PAGE 123 if (ibc==3) recyc_dist = sqrt( sq(x[ii]x[jj]) + sq(y[ii]y[jj]) + sq(z[ii]z[jj]) ); if (recyc_dist <= 0.5*(diam[ii]+diam[jj]) )// TRY THIS WITHOUT THE FOLLOWING:  y[i]<=1/tanhhar 1.0*(x[i] 0.5*hod) + 0.5*diam[i]/sinhhar  y[i] <= 1/tanhhar 1.0*(1.0*x[i]0.5*hod) + 0.5*diam[i]/sinhhar ) { contact=1; break; } } } if (count3 >= 1500) { count3=0; y_loc += diameter; cout <<"\nincreaseing y_loc\n"; cout <<"Recently placed particle < PAGE 125 if (neighbors[0][i]!=0) { cout <<"particle "< PAGE 126 n[3] = z[ii] z[jj]; n_magn = sqrt( n[1]*n[1] + n[2]*n[2] + n[3]*n[3] ); n[1] = n[1] / n_magn; n[2] = n[2] / n_magn; n[3] = n[3] / n_magn; temp2 = cx[ii]; g[1] = cx[ii] cx[jj] 0.5*( (diam[ii]*wy[ii]+diam[jj]*wy[jj])*n[3] (diam[ii]*wz[ii]+diam[jj]*wz[jj])*n[2] ); } // close if "normal" collision else // the collisions across a boundary { if (ibc==2 && idim==3) // collisions across the periodic hopper boundary { if (isup == 0) { cout << "\nCollision across the zpe riodic boundary!! Coll btn "< PAGE 127 g[2] = cy[ii] cy[jj] 0.5*( 1.0*(diam[ii]*wx[ii]+diam[jj]*wx[jj])*n[3] + (diam[ii]*wz[ii]+diam[jj]*wz[jj])*n[1] ); g[3] = cz[ii] cz[jj] 0.5*( (diam[ii]*wx[ii]+diam[jj]*wx[jj])*n[2] (diam[ii]*wy[ii]+diam[jj]*wy[jj])*n[1] ); } // close if j!=0 (pp coll) delta = overlap; g_dot_n = g[1]*n[1] + g[2]*n[2] + g[3]*n[3]; s[1] = g[1] g_dot_n*n[1]; s[2] = g[2] g_dot_n*n[2]; s[3] = g[3] g_dot_n*n[3]; s_magn = sqrt( sq(s[1]) + sq(s[2]) + sq(s[3]) ); if (s_magn<0.0000000000001) { s_magn=1; cout <<" \n\n AM I really HERE??? 2948 \n\n"; } s[1] = s[1] / s_magn; s[2] = s[2] / s_magn; s[3] = s[3] / s_magn; /* for (b=1; b<=3; b++)// for pw collisions, there were some very small +/flucs in the s vector in the normal // direction (limitation of machine precision?? due to e. g. g[2] g_dot_n n[2] should = 0, but is instead a very // small number). Hopefully this fix will mitigate that problem. { if ( fabs(s[b]) < 0.0000000001 ) { s[b] = 0; cout << \n\n AM I really HERE ?? 395 \n\n"; } } */ g_dot_s = g[1]*s[1] + g[2]*s[2] + g[3]*s[3]; // warn with output to screen if normal or tangential vectors are not of unit length sm = sqrt( sq(s[1]) + sq(s[2]) + sq(s[3]) ); nm = sqrt( sq(n[1]) + sq(n[2]) + sq(n[3]) ); if ( fabs(nm1) > 0.001 ) cout <<"\n\n\n\n\n WARNING!!! \n\n normal vector is not of unit length: \n\n "< PAGE 128 if (isup==0 ) { cout << setprecision(4) << setw(5) << "\nn = "< PAGE 129 //cout<<"diami is < PAGE 130 { SeparationDistance = (cur_sep(diam[ii]*0.5)(diam[jj]*0.5)); if (SeparationDis tance<=asperity) { SeparationDistance = asperity; } relcx=cx[ii]cx[jj]; relcy=cy[ii]cy[jj]; relcz=cz[ii]cz[jj]; RelativeVelocity=sqrt(((pow(relcx,2)+ (pow(relcy,2))+ (pow(relcz,2))))); hm_radius=((0.5*diam[ii])*(0.5*diam[ jj]))/((0.5*diam[ii])+(0.5*diam[jj])); Fn_Viscous_Magn =(6*3.14*viscosity_w*RelativeVelocity*hm_radius)*((hm_radius)/SeparationDistance); for (b=1; b<=3; b++) { Fnvfij[b] = Fn_Viscous_Magn n[b]; } } //Anshue void SoftPPCollisionModel(int ii, int jj) { // Define appropriate contact parameters for collision of interest: if (ii<=nLG) { if (jj<=nLG) { mu = mu_ll; k1 = k1_ll; k2 = k2_ll; kt = kt_f k1_ll; dd = dia; } else { mu = mu_sl; k1 = k1_sl; k2 = k2_sl; kt = kt_f k1_sl; dd = 0.5*(dia+dia_sm); } } if (ii>nLG) { if (jj<=nLG) { mu = mu_sl; k1 = k1_sl; k2 = k2_sl; kt = kt_f k1_sl; dd = dia_sm; cout <<"this should never occur, right?"; } else { mu = mu_ss; k1 = k1_ss; k2 = k2_ss; kt = kt_f k1_ss; dd = dia_sm; } } pp_total_integrations++; k_ratio = k2 / (k2k1); if ( g_dot_n < 0.0 ) // Loading { Fn_magn = k1 delta; } if ( g_dot_n <= 0.0 ) { delta_o = max( delta_o, 1.0 delta ( 1.0 k1 / k2 ) ); } ovlp_for_se = ovlp_for_se + sq(delta); if (delta>maxdelta) maxdelta=delta; if ( ( delta < delta_o && g_dot_n > 0.0 ) ) 130 PAGE 131 { part_collisions++; pp_collisions[ii]++; //collision counters pp_collisions[jj]++; pp_collisions_segn_stats++; for (k=0; k<=num_bin_o; k++) // Bin the overlap for p/o of distribution { if ( ( k_ratio delta_o / 0.5 / (diam[ii] + diam[jj]) 100 ) < bin_o[k][1] ) { // bin*[k][2] == count of velocities in bin k bin_o[k][2]++; break; } } sum_overlap = sum_overlap + k_ratio delta_o ; if(k_ratio delta_o < min_overlap) min_overlap = k_ratio delta_o; if(k_ratio delta_o > max_overlap) { // cout < PAGE 132 if (delta>0.1*diam[ii]) cout <<"\n X t "< PAGE 134 // cout << deltatorque= "<<0.015 Fn_magn wx[ii] / wi_magn<<" "<<0.015 Fn_magn wy[ii] / wi_magn<<" "<<0.015 Fn_magn wz[ii] / wi_magn< PAGE 135 torqueij_b[2]= Tib[ii][1][0]*torqueij_i[1] + Ti b[ii][1][1]*torqueij_i[2] + Tib[ii][1][2]*torqueij_i[3] ; torqueij_b[3]= Tib[ii][2][0]*torqueij_i[1] + Ti b[ii][2][1]*torqueij_i[2] + Tib[ii][2][2]*torqueij_i[3] ; torqueji_b[1]= Tib[jj][0][0]*torqueji_i[1] + Tib[jj][0][1]*torqueji_i[2] + Tib[jj][0][2]*torqueji_i[3] ; torqueji_b[2]= Tib[jj][1][0]*torqueji_i[1] + Tib[jj][1][1]*torqueji_i[2] + Tib[jj][1][2]*torqueji_i[3] ; torqueji_b[3]= Tib[jj][2][0]*torqueji_i[1] + Tib[jj][2][1]*torqueji_i[2] + Tib[jj][2][2]*torqueji_i[3] ; } else // transformation not required in 2D, but need to switch variable names anyway pw coll in xy plane looks fine, but in yz plane the spin has the wrong sign */ // { torqueij_b[1] = torqueij_i[1]; torqueij_b[2] = torqueij_i[2]; torqueij_b[3] = torqueij_i[3]; torqueji_b[1] = torqueji_i[1]; torqueji_b[2] = torqueji_i[2]; torqueji_b[3] = torqueji_i[3]; // } if(isup==0 ){ cout << setprecision(4)<<" \ntorqueij_b[1]b2="< PAGE 136 if (x[i] <= 0.5 (_hdiadiam[i]) && y[i]>0+y_origin) { wall_num=8; lbw[i]=wall_num; //NoLiquidBridges[i]=NoLiquidBridges[i]+1; ParticleWallLiquidBridge Model(i,wall_num,_hdia,_hod,_hhar); WallForceLiqBridge(i); } // Right Bin Wall if (x[i] >= 0.5 (_hdiadiam[i]) && y[i]>0+y_origin ) { wall_num=0; lbw[i]=wall_num; //NoLiquidBridges[i]=NoLiquidBridges[i]+1; ParticleWallLiquidBridge Model(i,wall_num,_hdia,_hod,_hhar); WallForceLiqBridge(i); } //Right hopper wall if ( ( _hha<90 && y[i]<1/_tanhhar 1.0*(x[i] 0.5*_hod) + 0.5*diam[i]/_sinhhar + y_origin && y[i] > 1.0*_tanhhar*(1.0*x[i] + 0.5*_hod) && y[i]>0+y_origin )  ( _hha==90.0 && y[i]<0.5*diam[i] && y[i]>y_origin && x[i]>0.5*_hod) ) { wall_num=1; lbw[i]=wall_num; //NoLiquidBridges[i]=NoLiquidBridges[i]+1; ParticleWallLiquidBridgeModel(i,wall_num,_hdia,_hod,_hhar); WallForceLiqBridge(i); } // Left Hopper Wall if ( (_hha<90 && y[i] < 1.0/_t anhhar (1.0*x[i]0.5*_hod ) + 0.5*diam[i]/_sinhhar + y_origin && y[i] > 1.0*_tanhhar*(x[i]+0.5*_hod) && y[i]>0+y_origin)  (_hh a==90.0 && y[i]<0.5*diam[i] && y[i]>y_origin && x[i] < 0.5*_hod) ) { wall_num=7; lbw[i]=wall_num; //NoLiquidBridges[i]=NoLiquidBridges[i]+1; ParticleWallLiquidBridge Model(i,wall_num,_hdia,_hod,_hhar); WallForceLiqBridge(i); } // right orifice wall //if (y[i]<=0 && x[i]>=0.5*(_hodd iam[i]) && y[i]>=0.17*_ho d0.025)// 'upper' hopper need not be included { //wall_num=3; //lbw[i]=wall_num; //NoLiquidBridges[i]=NoLiquidBridges[i]+1; 136 PAGE 137 // cout << "There has just been a collision bt n particle "<=0. 17*_hod0.025)// 'upper' hopper need not be included { //wall_num = 5; //lbw[i]=wall_num; //NoLiquidBridges[i]=NoLiquidBridges[i]+1; //cout << "There has just been a collision btn particle "< PAGE 138 if (SeparationDistance_wall<= 0.0) { SeparationDistance_wall=0.0; } Dimless_sep_wall = Sepa rationDistance_wall/(diam[ii]*0.5); NoLiquidBridges[ii]=NoLiquidBridges[ii]+1; //cout<<"Particle "< PAGE 139 Fn_Liq_Magn_wall = 3.14*diam[ii]*0.5*gamma_w*(exp((mikami_a_wall* Dimless_sep_wall)+mikami_b_wall)+mikami_c_wall); for (b=1; b<=3; b++) { Fnlbwij[b] = 1.0*Fn_Liq_Magn_wall n[b]; } } else { lbw[ii]=11; // particle is freed from liquid bridges NoLiquidBridges[ii]=NoLiquidBridges[ii]1; } } if (wall==1) { n[1] = 1.0 _coshhar; n[2] = _sinhhar; line= (_tanhhar*y[ii])x[ii]+(0.5*_hod); if (line < 0) line=line; SeparationDistance_wall = _coshhar*line (0.5*diam[ii]); if (SeparationDistance_wall<= 0.0) { SeparationDistance_wall=0.0; } Dimless_sep_wall = Sepa rationDistance_wall/(diam[ii]*0.5); NoLiquidBridges[ii]=NoLiquidBridges[ii]+1; //cout<<"Particle "< PAGE 140 n[1] = _coshhar; n[2] = _sinhhar; line= (_tanhhar*y[ii])+x[ii]+(0.5*_hod); if (line < 0) line=line; SeparationDistance_wall = _coshhar*line (0.5*diam[ii]); if (SeparationDistance_wall<= 0.0) { SeparationDistance_wall=0.0; } Dimless_sep_wall = Sepa rationDistance_wall/(diam[ii]*0.5); NoLiquidBridges[ii]=NoLiquidBridges[ii]+1; //cout<<"Particle "< PAGE 141 { mikami_a_wall = 1*1.9*(pow(Dimless_LiqBridge_Volume,0.51)); mikami_b_wall = (((1*0.016*(lo g(Dimless_LiqBridge_Volume )))0.76)*(pow(theta_wg,2))0.12*(log(Dimless_LiqBridge_Volume))+1.2); mikami_c_wall =(0.013*(log(Dimless_LiqBridge_Volume))) + 0.18; Fn_Liq_Magn_wall = 3.14*diam[ii]*0.5*gamma_w*(exp((mikami_a_wall* Dimless_sep_wall)+mikami_b_wall)+mikami_c_wall); for (b=1; b<=3; b++) { Fnlbwij[b] = 1.0*Fn_Liq_Magn_wall n[b]; } } else { lbw[ii]=11; // particle is freed from liquid bridges NoLiquidBridges[ii]=NoLiquidBridges[ii]1; } } if (wall==5) { n[1] = 1.0 ; n[2] = 0.0; SeparationDistance_wall = (x[ii] (0.5*_hod)(0.5*diam[ii])); if (SeparationDistance_wall<= 0.0) { SeparationDistance_wall=0.0; } Dimless_sep_wall = Sepa rationDistance_wall/(diam[ii]*0.5); NoLiquidBridges[ii]=NoLiquidBridges[ii]+1; //cout<<"Particle "< PAGE 142 } //Anshue void SoftPWCollisionModel(int ii, int wall) { // Define appropriate contact parameters for collision of interest: // if (wall==4) cout <<" in softpwcollmodel for wall=4 "; if (ii<=nLG) { muw = mu_lw; k1w = k1_lw; k2w = k2_lw; ktw = ktw_f k1_lw; dd = dia;} else { muw = mu_sw; k1w = k1_sw; k2w = k2_sw; ktw = ktw_f k1_sw; dd = dia_sm;} double deltao_wall = pwdeltao[wall][ii]; double pw_tan_dispx = pwtan_dispx[wall][ii]; double pw_tan_dispy = pwtan_dispy[wall][ii]; double pw_tan_dispz = pwtan_dispz[wall][ii]; pw_total_integrations++; if ( g_dot_n < 0.0 ) // Loading { Fn_magn = k1w delta; } if ( g_dot_n <= 0.0 ) { deltao_wall = max( deltao_wall, 1.0 delta ( 1.0 k1w / k2w ) ); } // if (wall==4) cout <<"\n n1, n2, n3, Fn_magn, delta deltao_wall = "< PAGE 143 if(k2w / (k2wk1w) deltao_wall < min_pwoverlap) min_pwoverlap = k2w / (k2wk1w) deltao_wall; if(k2w / (k2wk1w) deltao_wall > max_pwoverlap) { ///( cout < PAGE 144 psi_x = psi_total s[1]; psi_y = psi_total s[2]; psi_z = psi_total s[3]; psi_magn = sqrt( sq(psi_x) + sq(psi_y) + sq(psi_z) ); Ftij[1] = 1.0 psi_x / psi_magn min( ktw psi_magn, muw Fn_magn); Ftij[2] = 1.0 psi_y / psi_magn min( ktw psi_magn, muw Fn_magn); Ftij[3] = 1.0 psi_z / psi_magn min( ktw psi_magn, muw Fn_magn); for (b=1; b<=3; b++) Ftji[b] = 1.0 Ftij[b]; if ((k1w < 1500 && delta > 0.5*dd)  (k1w < 15000 && k1w > 1500 && delta > 0.3*dd)  (k1w > 15000 && delta > 0.2*dd)  (k1w > 15000 && k1w < 150000 && delta > 0.15*dd)  ( k1w > 150000 && delta > 0.1*dd)) // && t>=topen) { overlap_warning++; if (countlargedelta < 100) { countlargedelta++; cout << setw(7)< PAGE 145 torqueij_i[1] = torqueij_i[1] mu_rw Fn_magn wx[ii] / wi_magn; torqueij_i[2] = torqueij_i[2] mu_rw Fn_magn wy[ii] / wi_magn; torqueij_i[3] = torqueij_i[3] mu_rw Fn_magn wz[ii] / wi_magn; } } if(isup==0 ){ cout << setprecision(4)<<"\nFnij[1] = "< PAGE 146 torqueji_b[1] = torqueji_i[1]; torqueji_b[2] = torqueji_i[2]; torqueji_b[3] = torqueji_i[3]; // } if(isup==0 ){ cout << setprecision(4)<<" \ntorqueij_b[1]b2="< PAGE 147 */ ifstream input("seg8r_ input.txt", ios::in); input >> np; input.ignore(80, '/'); input >> ifill; input.ignore(180, '/'); input >> randit; input.ignore(160, '/'); input >> nreps; input.ignore(160, '/'); if (nreps>100) {cout <<"\nSorry, the input 'nreps cannot exceed 100 replica tions\n"; return(57);} input >> ibc; input.ignore(160, '/'); if (ibc==0) {cout <<"\nSorry, the option 'periodic shear flow' is not supported in this version of the code\n"; return(22);} input >> nrecyc; input.ignore(160, '/'); input >> irot; input.ignore(120, '/'); input >> isearch; input.ignore(120, '/'); input >> idim; input.ignore(120, '/'); input >> z_width; input.ignore(80, '/'); input >> uho; input.ignore(80, '/'); input >> uhht; input.ignore(80, '/'); input >> uhdia; input.ignore(80, '/'); input >> uhod; input.ignore(80, '/'); input >> uhha; input.ignore(80, '/'); input >> hht; input.ignore(80, '/'); input >> hdia; input.ignore(80, '/'); input >> hod; input.ignore(80, '/'); input >> hha; input.ignore(80, '/'); input >> idiam; input.ignore(90, '/'); 147 PAGE 148 input >> idd; input.ignore(160, '/'); input >> dia; input.ignore(160, '/'); input >> sd; input.ignore(160, '/'); input >> dens; input.ignore(90, '/'); input >> nLG; input.ignore(90, '/'); if (nLG>np && idiam == 2) { cout << "\n\n Error! Number of large particles (nLG) is greater than the total number of particles (np)!\n"; return(22);} input >> dia_sm; input.ignore(160, '/'); input >> sd_sm; input.ignore(160, '/'); input >> dens_sm; input.ignore(90, '/'); input >> num_sam; input.ignore(160, '/'); input >> sam_rad; input.ignore(160, '/'); input >> dt; input.ignore(80, '/'); input >> itadj; input.ignore(80, '/'); input >> tfinal; input.ignore(80, '/'); input >> twrite; input.ignore(120, '/'); input >> movname; input.ignore(120, '/'); input >> mu_ll; input.ignore(80, '/'); input >> mu_ss; input.ignore(80, '/'); input >> mu_sl; input.ignore(80, '/'); input >> mu_r; input.ignore(80, '/'); input >> mu_lw; input.ignore(80, '/'); input >> mu_sw; input.ignore(80, '/'); input >> mu_rw; input.ignore(80, '/'); input >> k1_ll; input.ignore(80, '/'); 148 PAGE 149 input >> k2_ll; input.ignore(80, '/'); input >> k1_ss; input.ignore(80, '/'); input >> k2_ss; input.ignore(80, '/'); input >> k1_sl; input.ignore(80, '/'); input >> k2_sl; input.ignore(80, '/'); input >> kt_f; input.ignore(80, '/'); input >> k1_lw; input.ignore(80, '/'); input >> k2_lw; input.ignore(80, '/'); input >> k1_sw; input.ignore(80, '/'); input >> k2_sw; input.ignore(80, '/'); input >> ktw_f; input.ignore(80, '/'); input >> ivdamp; input.ignore(80, '/'); input >> igrav; input.ignore(80, '/'); input >> gravacc; input.ignore(80, '/'); input >> ivdist; input.ignore(80, '/'); input >> ifile; input.ignore(80, '/'); input >> imovie; input.ignore(80, '/'); input >> ihe; input.ignore(80, '/'); input >> isup; input.ignore(80, '/'); input >> ipos; input.ignore(80, '/'); input >> icoll; input.ignore(80, '/'); input >>gamma_w; input.ignore(80,'/'); input >>viscosity_w; input.ignore(80,'/'); input >>theta_wg; input.ignore(80,'/'); input >>theta_ww; input.ignore(80,'/'); input >>Volume1; 149 PAGE 150 input.ignore(80,'/'); input >>Volume2; input.ignore(80,'/'); input >>asperity; input.ignore(80,'/'); input >>Volume_Wall; input.ignore(80,'/'); input.close(); /* ifstream get_exec("ss", ios::in); get_exec.ignore(80, '.'); get_exec.ignore(80, '.'); get_exec>>exec_name; get_exec.close(); */ ifstream get_exec("ss", ios::in); get_exec.ignore(280, '.'); get_exec.ignore(280, '.'); get_exec.ignore(280, '.'); get_exec>>exec_name; get_exec.close(); cout << "Executable name: \t."< PAGE 151 //Anshus NoLiquidBridges = new int[np+1]; for (i=1;i<=np;i++) NoLiquidBridges[i]=0; lbw = new int[np+1]; for(i=1;i<=np;i++) lbw[i]=11; //Anshue x = new double[np+1]; y = new double[np+1]; z = new double[np+1]; cx = new double[np+1]; cy = new double[np+1]; cz = new double[np+1]; wx = new double[np+1]; wy = new double[np+1]; wz = new double[np+1]; mass = new double[np+1]; diam = new double[np+1]; moment = new double[np+1]; eta = new double[np+1]; xi = new double[np+1]; lambda = new double[np+1]; w_magn = new double[np+1]; ave_velocity = new double[np+1]; pw_nn_list = new int[np+1]; skip_part = new int[np+1]; skip_search = new int[np+1]; pp_collisions = new int[np+1]; pw_collisions = new int[np+1]; discharget = new int[np+1]; Fnsum = new double*[4]; for (i=0; i<=3; i++) Fnsum[i] = new double[np+1]; torquesum = new double*[4]; for (i=0; i<=3; i++) torquesum[i] = new double[np+1]; discharge = new double*[9]; for (i=0; i<=8; i++) discharge[i] = new double[np+1]; neighbors = new int*[700]; for (i=0; i<=699; i++) neighbors[i] = new int[np+1]; pwdeltao = new double*[9]; for (i=0; i<=8; i++) pwdeltao[i] = new double[np+1]; pwtan_dispx = new double*[9]; for (i=0; i<=8; i++) pwtan_dispx[i] = new double[np+1]; pwtan_dispy = new double*[9]; for (i=0; i<=8; i++) pwtan_dispy[i] = new double[np+1]; pwtan_dispz = new double*[9]; 151 PAGE 152 for (i=0; i<=8; i++) pwtan_dispz[i] = new double[np+1]; /*RUN MULTIPLE RANDOM ITERATIONS (using do...while loop) */ do{ time_t tm1 = time(NULL); // note the st art time for the file "*_report_out.dat" time_t tm10 = time(NULL); // note the start time for the periodic updates of the system status time_t tm20; double del_time; cout <<"Start time stamp: \t"< PAGE 153 cout << setw(10) << setprecision(4); cout << "itadj << itadj < PAGE 154 } if (ibc == 3) { pw_max_dist = 1000.0; for (i=1; i<=np; i++) pw_nn_list[i] = 0; } dtinv = 1.0 / dt; npinv = 1.0 / np; if ( mu_r != 0  mu_rw !=0 ) irf = 1; else irf = 0; grav = 0.0; if (igrav == 1) grav = gravacc; filename1 = movname + "_hdp" + itos(curr_randit) + "_out.dat"; filename2 = movname + "_time" + itos(curr_randit) + "_distn_out.dat"; filename3 = movname + "_flow" + itos(curr_randit) + "_out.dat"; filename4 = movname + "_startup_pos" + itos(cu rr_randit) + ".dat"; filename33 = movname + "_startup_cont" + itos(curr_randit) + ".dat"; filename5 = movname + "_mass" + itos(curr_randit) + "_distn_out.dat"; filename6 = movname + "_par_states.dat"; filename7 = movname + "_spheresandhopper.xml"; filename8 = movname + "_standard_colormap.dat"; filename9 = movname + "_RTD_colormap.dat"; filename10 = movname + "_fill_part_states.dat"; filename12 = movname + "_temp_pos.dat"; filename32 = movname + "_temp_contacts.dat"; filename13 = movname + "_initfill_part_states.dat"; filename16 = movname + "_H25_velocity_distn.dat"; filename17 = movname + "_H12_velocity_distn.dat"; filename18 = movname + "_H08_velocity_distn.dat"; filename21 = movname + "_striped_colormap.dat"; filename22 = movname + "_pw_nn_colormap.dat"; filename23 = movname + "_report_out.dat"; filename24 = movname + "_summary_massdistn_out.dat"; filename25 = movname + "_ave_velocity_colormap.dat"; filename26 = movname + "_sifting_colormap.dat"; filename27 = movname + "_sifting_velocity_vectors.dat"; // filename28 = movname + "_local_xi_colormap.dat"; filename29 = movname + "_tecplot_time_mass_frac.dat"; filename30 = movname + "_fixed_RTD_colormap.dat"; filename31 = movname + "_diamass.dat"; filename34 = movname + "_liquidbridge_colormap.dat"; filename35 = movname + "_mydata1.dat"; filename36 = movname + "_mydata2.dat"; hhar = hha Pi / 180; uhhar = uhha Pi / 180; 154 PAGE 155 coshhar = cos(hhar); sinhhar = sin(hhar); tanhhar = tan(hhar); tanuhhar = tan(uhhar); cout<<"\nhopper angles and trig fcns (cos, tan, sin): "< PAGE 156 qq = HRand(.1); rr = sq(pp) + sq(qq); } while (rr >= 1.0); rr = sqrt(2.0 log(rr) / rr); if (idd == 2) diam[i] = 1.0 sd pp rr + dia; if (idd == 3) diam[i] = exp(1.0 sd_ln pp rr + dia_ln); } while (diam[i] < 0.02); } pp_collisions[i] = 0; /* tempcx[i] = 0.0; tempcy[i] = 0.0; tempcz[i] = 0.0; oldtempcx[i] = 0.0; oldtempcy[i] = 0.0; oldtempcz[i] = 0.0; */ for (j=1; j<=np; j++) { cur_sep = 10.0; } } if (idd==1 && idiam==2 ) { for (i=1; i<=nLG; i++) { diam[i] = dia; } for (i=nLG+1; i<=np; i++) { diam[i] = dia_sm; } mindiam=dia_sm; maxdiam=max(dia,dia_sm); } if (idd!=1 && idiam==1) { for (i=1; i<=np; i++) { avediam = avediam + diam[i]; } avediam = avediam / np; } 156 PAGE 157 if (idd!=1 && idiam==2 ) // binary, gaussian distribution { for (i=1; i<=nLG; i++) { do { do { pp = HRand(.1); qq = HRand(.1); rr = sq(pp) + sq(qq); } while (rr >= 1.0); rr = sqrt(2.0 log(rr) / rr); if (idd==2) diam[i] = 1.0 sd pp rr + dia; if (idd==3) diam[i] = exp(1.0 sd_ln pp rr + dia_ln); } while (diam[i] < 0.02); } for (i=nLG+1; i<=np; i++) { do { do { pp = HRand(.1); qq = HRand(.1); rr = sq(pp) + sq(qq); } while (rr >= 1.0); rr = sqrt(2.0 log(rr) / rr); if (idd==2) diam[i] = 1.0 sd_sm pp rr + dia_sm; if (idd==3) diam[i] = exp( 1.0 sd_sm_ln pp rr + dia_sm_ln); } while (diam[i] < 0.02); } } if (idd!=1 && idiam==2) // calculate averag e diameters for binary, gaussian distribution { for (i=1; i<=nLG; i++) { avediamLG = avediamLG + diam[i]; maxdiam = max(maxdiam, diam[i]); } for (i=nLG+1; i<=np; i++) 157 PAGE 158 { avediamSM = avediamSM + diam[i]; mindiam = min(mindiam, diam[i]); } avediamSM = avediamSM / (npnLG); avediamLG = avediamLG / nLG; } if (diam[i] PAGE 159 } if (idim==2) { minmass = dens_sm*dia_sm*dia_sm*Pi*0.25; maxmass = dens*dia*dia*Pi*0.25; } } if (idd!=1) { if (idim==3) { minmass = dens_sm*avediamSM*avediamSM*avediamSM*Pi*0.166666667; maxmass = dens*avediamLG*avediamLG*avediamLG*Pi*0.166666667; } if (idim==2) { minmass = dens_sm*aved iamSM*avediamSM*Pi*0.25; maxmass = dens*avediamLG*avediamLG*Pi*0.25; } } x_fines = minmass*(npnLG)/(minmass*(npnLG) + maxmass*nLG); } loopstop = min(50, nLG); loopstop2 = min(50, np); if (np<50) loopstop=loopstop2=np; if (idiam==1) { for (i=1; i<=loopstop2; i++) { cout < PAGE 160 for (i=nLG+1; i<=nLG+1+loopstop2; i++) // small particles { cout < PAGE 161 input.ignore(1); input >> cx[i]; input.ignore(1); input >> cy[i]; input.ignore(1); input >> cz[i]; input.ignore(1); if (i==1  i%100==0) cout <<"part. "<> i; input_cont.ignore(1); input_cont >> j; input_cont.ignore(1); input_cont >> cur_sep; input_cont.ignore(1); input_cont >> delta_o; input_cont.ignore(1); input_cont >> tan_dispx; input_cont.ignore(1); input_cont >> tan_dispy; input_cont.ignore(1); input_cont >> tan_dispz; input_cont.ignore(1); cout <cur_sep = cur_sep; a>delta_o = delta_o; a>tandx = tan_dispx; a>tandy = tan_dispy; a>tandz = tan_dispz; } while (input_cont.peek()!=1); if (irot==1  irf==1) { for (i=1; i<=np; i++) { wx[i] = 0.25 ( HRand(rand_factor) ); wy[i] = 0.25 ( HRand(rand_factor) ); 161 PAGE 162 wz[i] = 0.25 ( HRand(rand_factor) ); } } for (i=1; i<=np; i++) // account for particles starting 'below' hopper { if (y[i] < 0.25*hod) { skip_search[i]=1; topen = 0; } if (y[i] < 2.0*hod) skip_part[i]=1; if (ibc==3) { pw_nn_list[i]=0; xz_sq = sqrt( sq(x[i]) + sq(z[i]) ); if ( ( y[i] > 1.0/tanhhar*(0.5*hod+0.5*hdia2.500*dia) + 2.500*dia/sinhhar && xz_sq > 0.5*hdia 2.500*dia )  ( y[i] < 1.0/tanhhar*(0.5*hod+0 .5*hdia2.500*dia) + 2.500*dia/sinhhar && y[i] < 1.0/tanhhar*xz_sq 0.5*hod/tanhhar + 2.500*dia/sinhhar )  y[i] < 2.500*dia ) pw_nn_list[i]=1; // particle is NOT in the near neighbor zone i.e. it is in the central region } } } // close if ifill==4 if (ifill==5) //Random HCP fill pl aces particles on a hexagonal close p acked lattice, then randomly shift particles from // these lattice points. Useful for systems with large solids fractions. Only to be used in 2D shear. { spacing = ceil( 1.0 np / floor( sqrt(1.0*np)) ); spacing = 2.0 / sqrt(3.0) ymax / (spacing+1.0); x[1] = 0.0; y[1] = 0.0; xinc = 0.5 spacing; yinc = 0.5 sqrt(3.0) spacing; old_yinc = 1.0 yinc; for (i=2; i<=np; i++) { x[i] = x[i1] + xinc; y[i] = y[i1] + yinc; if ( (x[i] >= ymax spacing && yinc < 0.0)  (x[i] >= ymax 0.5 spacing && yinc > 0.0) ) { if (yinc < 0.0) { x[i] = 0.0; y[i] = y[i1] yinc; 162 PAGE 163 } else { x[i] = 0.0; y[i] = y[i1] + 2.0 yinc; aa = old_yinc; old_yinc = yinc; yinc = aa; } } if (y[i] > 5.0) { x[i] = x[i] + xinc; y[i] = y[i] yinc; do { i++; x[i] = x[i1] + 2.0 xinc; y[i] = y[i1]; } while (i<=np); } aa = old_yinc; old_yinc = yinc; yinc = aa; } HRand(.1); HRand(.1); HRand(.1); HRand(.1); HRand(.1); for (i=1; i<=np; i++) //add random adjustment to position { x[i] = x[i] + HRand(.1) 0.5 (spacing diam[i]); y[i] = y[i] + HRand(.1) 0.5 (spacing diam[i]); } } // THE "NEW" RANDOM FILL ================================================================================== ========== if (ifill==1) { if (idiam==1  idiam==2) { if (ibc==2ibc==3) RandomFill(0, hht, hdia, hod, hha); 163 PAGE 164 if (ibc==4) RandomFill_Pharma(0, hht, hdia, hod, hha); } // close if idiam == 1 or 2 } // close if ifill=1 if (ifill==9) // trickle fill from hopper located above 'main' hopper now using ASTM tester geometry { // Fill upper hopper: RandomFill(uho, uhht, uhdia, uhod, uhha); for (i=1; i<=np; i++) { cx[i] = 0.5*HRand(0.1); cy[i] = 1.0 fabs(0.8*HRand(0.1)); cz[i] = 0.15*HRand(0.1); } for (i=1; i<=np; i++) { pw_nn_list[i] = 1; } // d on't use the pw_nn_list for ifill==9...could be used but would need to modify the code then if (irot==1  irf==1 ) { for (i=1 ; i<=np; i++) {wx[i] = 0.1*HRand(0.1 ); wy[i] = 0.1*HRand(0.1); wz[i] = 0.1*HRand(0.1); }} t = 0; search_count[2] = 10; // ensures NN search at t=1 dt = 1.5*dt; //artificially increase time step size for the initial filling of upper hopper ONLY! // Start mini time loop: do { t++; // PP Collision search if (idiam==1 && search_count[2]>=0.55*dia) { for (i=1; i<=np; i++) { neighbors[0][i] = 0;} // reset the number of neighbors column to zero maxneighbors=0; ave_neighbors = 0; FindNeighbors( 1.55*dia, 1, np, 1, np); //search over all particles, for all possible neighbors search_count[2] = 0; // PrintNeighborList(); } if (idiam==2 && (search_count[1]>0.55*dia  search_count[2]>0.5*0.55*(dia+dia_sm)  search_count[4]>0.55*dia_sm)) { //cout << "search count = "< PAGE 168 for (i=1; i<=np; i++) { temp_out< PAGE 169 z[1] = 0.0; cx[1] = 0.4; cy[1] = 0.0; cz[1] = 0.0; wx[1] = 0.0; wy[1] = 0.0; wz[1] = 50.0; */ /* x[2] = 2; y[2] = 2; z[2] = 0; cx[2] = 0; cy[2] = 534134; cz[2] = 0; */ //plot of PSI1 and PSI2 on flat bottomed hopper wall /* x[1] = 1.0; // xy plane y[1] = .051; z[1] = 0.0; cx[1] = 1.0*sin(gravacc*Pi/180); cy[1] = 1.0*cos(gravacc*Pi/180); cz[1] = 0.0; */ /* z[1] = 0; // yz plane y[1] = .051; x[1] = 1.0; cz[1] = 1.0*sin(gravacc*Pi/180); cy[1] = 1.0*cos(gravacc*Pi/180); cx[1] = 0.0; */ x[1] = 1; y[1] = .6; z[1] = 0; cx[1]=8.0; cy[1] = 0; cz[1] = 0.001; x[2] = 3.; y[2] = .6; z[2] = 0.0; cx[2] = 8; cy[2] = 0.001; cz[2] = 0; if (idim==2) { z[1] = z[2] = 0.0; cz[1] = cz[2] = 0.0; wx[1] = wx[2] = 0.0; wy[1] = wy[2] = 0.0; } 169 PAGE 170 if (idim==3) { wx[1] = wx[2] = 0.0; wy[1] = wy[2] = 0.0; wz[1] = wz[2] = 0.0; } } /*ASSIGN RANDOM VELOCITY VECTORS TO EACH PARTICLE */ if (ifill == 1  ifill == 3  ifill == 5  ifill == 7  ifill == 8 ) // i.e. Don't run if fill is hard coded or hopper fil l { for (i=1; i<=np; i++) { cx[i] = 0.5 ( HRand(0.1) ); cy[i] = 1.0 fabs( 0.8 ( HRand(0.1)) ); cz[i] = 0.15 ( HRand(0.1) ); } if (irot==1  irf==1 ) { for (i=1; i<=np; i++) { wx[i] = 0.25 ( HRand(rand_factor) ); wy[i] = 0.25 ( HRand(rand_factor) ); wz[i] = 0.25 ( HRand(rand_factor) ); } } } if (idim == 2) // if 2_d reset z velocities to zero { for (i=1; i<=np; i++) { cz[i] = wx[i] = wy[i] = 0.0; maxzvel=0.0; } } for (i=1; i<=loopstop; i++) { cout << setw(6)< PAGE 171 { cout <<"\nx_fines = "< PAGE 172 sandh_out << "\t PAGE 173 sandh_out << "\t\t\t PAGE 174 if (ifill==9) { sandh_out << "\t\t PAGE 175 sandh_out << "\t\t PAGE 176 sandh_out << "\t\t\t PAGE 177 alph = 0.5*(hdiahod); for (b = 0; b<=71; b++) /// was <=35 { sinb = sin(2.50*b*Pi/180.0); /// was 5.0 instead of 2.5 in each of these four rows cosb = cos(2.50*b*Pi/180.0); sinbp1 = sin(2.50*(b+1)*Pi/180.0); cosbp1 = cos(2.50*(b+1)*Pi/180.0); // instead of 35.5 it was 0.5*(hdiahod)/tanhhar //cylinder sandh_out << "\t\t PAGE 178 }// close if ibc==4 sandh_out << "\t PAGE 179 num_bin_x_sift = round(hdia/3/dia); // the 3 > approx 3 large particle diameters num_bin_y_sift = round(2*hht/3/dia); deltabin_sift = hdia / num_bin_x_sift; //cout << " \nnum_bin_x_sift="< PAGE 180 psd_out << setprecision(5) << bin_psd[k][1] <<" "<<0.5*(bin_psd[k][1]+bin_psd[k1][1])<<" "<< setprecision(6) << bin_psd[k][2]<<" << bin_psd[k][2]/sum_of_psd/deltabin_psd< PAGE 181 } for (k=0; k<=num_bina; k++) { binax[k][1] = 2.5 + k deltabina; binay[k][1] = 2.5 + k deltabina; } /*CALCULATE SOLID AND VOID AREAS */ sarea = 0; if (ibc==2  ibc==3) // not rigorous... fix or delete { ymax = 0.0; if (idim==2) { for (i=1; i<=np; i++) { sarea = sarea + sq(diam[i])*Pi*0.25; ymax = max(ymax, y[i]); } totalarea = ymax hdia sq(0.5*(hdiahod))*tanhhar; varea = totalarea sarea; sfrac = sarea / totalarea; } if (idim == 3) { svol = 0.0; for (i=1; i<=np; i++) { svol = svol + 0.1666667*Pi*diam[i]*diam[i]*diam[i]; ymax = max(ymax, y[i]); } if (ibc==2) { totalvol = z_width*(ymax hdia sq( 0.5*(hdiahod) ) / tanhhar); } if (ibc==3) totalvol = Pi*sq(hdia/2)*(ymax 0.5*(hdiahod)/tanhhar) + Pi/3*(sq(hdia/2)*(hdia/2/tanhhar) sq(hod/2)*(hod/2/tanhhar)); vvol = totalvol svol; sfrac = svol / totalvol; } } cout <<" \ninitial solids fraction (approx!)="< PAGE 182 for (t=1; t<=tfinal; t++) { //cout<<"This is time Loop +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"< PAGE 183 cur_sep = min( sqrt( sq(x[i]x[j]) + sq(y[i]y[j]) + sq(z[i]z[j]) ), sqrt( sq(x[i]x[j]) + sq(y[i]y[j]) + sq(tempzz[j]) ) ) ; } if (cur_sep < 0.5*(diam[i] + diam[j])) // separtn[i][j] = min( separtn[i][j], cur_sep ); { // separtn[0][n_col] = i; // separtn[1][n_col] = j; // separtn[2][n_col] = cur_sep; n_col++; if (n_col >= 6*np 3 ) { cout <<" warning!! number of collisions has exceeded length of separation matrix"; return(1); } } // close if cur_sep } // closes for j } // closes skip_search }//closes for i } // closes if isearch == 1 /*NEARNEIGHBOR COLLISION SEARCH */ /*NEARNEIGHBOR COLLISION SEARCH PART A: DETERMINATION OF NEIGHBORS */ if (isearch==3) { // cout <<"\n entering if isearch==3 : "; if (is_recycled == 1) // initiates an u pdate of neighbor lists if a particle has been recycled { is_recycled = 0; search_count[2] = 10.0; } // cout <<"\nreseting numberof neighbours to zero : "; // argument list for FindNeighbors(radius, starting pa rticle index for particle to find neighbors of, final particle index for particle to find neighbors of, starting index of neighbors to look for, final index of neighbors to look for) if (idiam==1 && search_count[2]>=0.55*dia) 183 PAGE 184 { n_updates++; for (i=1; i<=np; i++) { neighbors[0][i] = 0; // reset the number of neighbors column to zero } maxneighbors=0; ave_neighbors = 0; FindNeighbors( 0.55*dia+dia, 1, np, 1, np); //search over all particles, for all possible neighbors search_count[2] = 0; // cout <<"\n Neighbors updated at t = "< PAGE 186 { //cout<<"index of Liquid Bridge map is ::"<time is :: "< PAGE 187 //cout<<"Volumej is < PAGE 188 for (i=1; i<=np; i++) { for (j=1; j<=np; j++) { // cout< PAGE 189 n_col = 0; if (t>topen) gate = 1; // slide gate is open else gate = 0; // slide gate is closed ParticleWallCollisions(0, hht, hdia, hod, hha, gate ); /*PARTICLE WALL LIQUID BRIDGES */ //Anshus: Call the Particle Wall Liquid Bridges Function fo r modifying with liquidbridge force with the wall. ParticleWallLiquidBridges(0,hht,hdia,hod,hha,gate); //Anshue //Anshus /*PARTICLEPARTI CLE LIQUID BRIDGES */ for(i=1;i<=np;i++) for(l=LiquidBridge_map[i];l;l=l>nex) { if (l>time==t) //cout<<"The code has entered the Particle Particle lIQUID bRIDGE "< PAGE 190 { // remove ?? the y>2.500*dia to prevent particles from sneaking through the slide gate during filling?? test_vel = sq(cx[i]) + sq(cy[i]) + sq(cz[i]); if ( test_vel > sq(pw_maxvel) ) pw_maxvel = sqrt( test_vel ); } } // calculate max possible distance traveled pw_max_dist += pw_maxvel dt; //cout < PAGE 191 // cout <<"\nParticle "<=topen) // hopper is 'open' { avedischtime++; // counts timesteps after discharge begins for the "by time" outflow distributions //Anshudischarges if (avedischtime%500 == 0) { // output avedischtime*dt for (i=1;i<=np;i++) { //cout< PAGE 192 if (islarge==0) { if (i<=nLG) { cout <<"\nt = "< PAGE 193 tavedischdiam = tavedischdiam + diam[i]; tavedischmass = tavedischmass + mass[i]; if (i<=nLG) tavedischmassL = tavedischmassL + mass[i]; if (i>nLG) tavedischmassS = tavedischmassS + mass[i]; if (idim == 2) tavedischdens = tavedischdens + 4*mass[i]/Pi/diam[i]/diam[i]; if (idim == 3) tavedischdens = tavedischdens + 6*mass[i]/Pi/diam[i]/diam[i]/diam[i]; for (k=0; k<=num_bin_hdp; k++) // bins the position of particle discharges { if ( x[i] < bin_hdp[k][1]) // bin*[k][1] == bin limit { // bin*[k][2] == count of positions in bin k bin_hdp[k][2]++; bin_hdp[k][5] += cy[i]; if ((i<=nLG && idiam==2)  idiam==1) // count large particles bin_hdp[k][3]++; if (i>nLG && idiam==2) // count small particles bin_hdp[k][4]++; break; } } } if ( ( np<500 && avedischtime==0.1/dt )  ( np>=500 && avedischtime == 10*twrite && ibc>=3 )  (np>=500 && avedischtime == 5*twrite && ibc==2 ) ) //used to be ==2000 { if (tdistnlocation > 198) cout <<"\n\nWARNING tdistnlocation > 200... impending seg fault!"< PAGE 194 tdistnlocation++; if (tdistnlocation > 198) cout <<"\n\nWARNING tdistnlocation > 200... impending seg fault!"< PAGE 195 } //close for i } //close if t>=topen /*WRAPUP FOR SOFTPARTICLE MODEL */ for (i=1; i<=np; i++) { if (skip_part[i] == 0) { //cout <<"\n i, mass[i], Fnsum1i, cx_pre, cx_post = < PAGE 196 Fnsum[b][i] = 0.0; torquesum[b][i] = 0.0; } } /*INTEGRATE FOR POSITION AND ANGULAR DISPLACEMENT VECTORS */ for (i=1; i<=np; i++) //later, worry about better integration schemes, etc...adjustable time, implicit, explicit T.I. { if (skip_part[i]==0) { x[i] = x[i] + dt cx[i]; y[i] = y[i] + dt cy[i]; z[i] = z[i] + dt cz[i]; } if (idim==3 ) { eta_dot = wx[i] cos(lambda[i]) / cos(xi[i]) + wy[i] 1.0 sin(lambda[i]) / cos(xi[i]); xi_dot = wx[i] sin(lambda[i]) + wy[i] cos(lambda[i]); lambda_dot = wx[i] 1.0 tan(xi[i]) cos(lambda[i]) + wy[i] tan(xi[i]) sin(lambda[i]) + 1.0 wz[i]; eta[i] = eta[i] + dt eta_dot; xi[i] = xi[i] + dt xi_dot; lambda[i] = lambda[i] + dt lambda_dot; if (eta[i]>2.0*Pi) eta[i] = eta[i] 2.0*Pi; if (xi[i]>2.0*Pi) xi[i] = xi[i] 2.0*Pi; if (lambda[i]>2.0*Pi) lambda[i] = lambda[i] 2.0*Pi; // cout < PAGE 197 } //close for i /*PERIODIC BOUNDARIES (FOR PLANE FLOW HOPPER) */ if (ibc==2) { for (i=1; i<=np; i++) { if (z[i] < 0.5*z_width  z[i] > 0.5*z_width) { corz=(int)(floor(2.0*z[i]/z_width+0.5)); z[i] = z[i] 1.0*corz*z_width; } } } /*CALCULATE CURRENT SYSTEM ENERGY */ pe = 0.0; ke = 0.0; re = 0.0; te = 0.0; ve = 0.0; se = 0.5 k1 ovlp_for_se; //if (se > 100000) cout <<"\nse calc: t se k1 ovlp_for_se "< PAGE 198 /*SET SKIP_SEARCH[i] & SKIP_PART[i], OR RECYCLE PARTICLE S, IF REQUIRED */ for (i=1; i<=np; i++) { if (y[i] < 1.0*hod && skip_search[i]==0) { skip_search[i]=1; //cout<<"Particle "<nLG) { cout <<"\nt = "< PAGE 199 /* surf_det_count = 0; contacts = 0; // if (recyc_events>=25) if (num_recycled<=1  num_recycled==25  num_recycled==75  num_recycled==200  num_recycled==350  num_recycled==500  num_recycled==800  num_r ecycled==1400  num_recycled==2500  num_recycled==3500  nu m_recycled==4900) { // recyc_events = 0; surface_min = 1.5*hod; // "reset" surface locations surface_max = 1.2*hht; do { do Rnp = 1 + floor(0.5*np*(HRand(.1)+1.0)); while (y[Rnp]<1.5*hod  z[Rnp]>0.4*z_width  z[Rnp]<0.4*z_width); // neglect particles in neighborhood of outer edge or the orif ice (fewer contacts, but still in "bulk") surf_det_count++; for (ir=1; ir<=Rnp1; ir++) // use 'ir' instead of 'i' (already in a for 'i' loop { if (collide[ir][Rnp]==1) contacts++; } for (j=Rnp+1; j<=np; j++) { if (collide[Rnp][j]==1) contacts++; } if (contacts > 1) surface_min = max(surface_min, y[Rnp]); if (contacts <= 1) surface_max = min(surface_max, y[Rnp]); } while (surf_det_count < floor(0.1*np));//e xamine 10% of particles to determine surface location recycle_ht = 0.5 (surface_max + surface_min) + max(1.0,surface_maxsurface_min); // last term MAY result in large y ht (?) cout <<"\n\nNEW recycle ht = "< PAGE 200 { j=j+5; test_dist = sqrt( sq(k*hdia/6.0 x[j]) + sq(y_test y[j]) ); // cout < PAGE 201 contact = 0; is_recycled=1; num_recycled++; discharget[i]=large_num; // reset discharge indicator cout <<"\n Particle "< maxvel*maxvel ) { maxvel = sqrt( cc ) ; fastest_p = i; } if ( fabs(cx[i]) > fabs(maxxvel)) { maxxvel = cx[i]; fastest_x_p=i; } if ( fabs(cy[i]) > fabs(maxyvel)) { maxyvel = cy[i]; fastest_y_p=i; } if ( fabs(cz[i]) > fabs(maxzvel)) { maxzvel = cz[i]; fastest_z_p=i; } if ( ss > maxrot*maxrot ) { maxrot = sqrt( ss ); fastest_rot_p=i; } } } } if (idiam==2) // find maximum velocities for large and sm all particles....for use in near neighbor searching { for (i=1; i<=nLG; i++) { if (skip_search[i]==0) { 201 PAGE 202 cc = sq(cx[i]) + sq(cy[i]) + sq(cz[i]); ss = sq(wx[i]) + sq(wy[i]) + sq(wz[i]); if ( cc > maxvel_L*maxvel_L ) { maxvel_L = sqrt( cc ) ; fastest_p_L = i; } if ( ss > maxrot*maxrot ) { maxrot = sqrt( ss ) ; fastest_rot_p = i; } } } for (i=nLG+1; i<=np; i++) { if (skip_search[i]==0) { cc = sq(cx[i]) + sq(cy[i]) + sq(cz[i]); ss = sq(wx[i]) + sq(wy[i]) + sq(wz[i]); if ( cc > maxvel_S*maxvel_S ) { maxvel_S = sqrt( cc ) ; fastest_p_S = i; } if ( ss > maxrot*maxrot ) { maxrot = sqrt( ss ) ; fastest_rot_p = i; } } } maxvel = max(maxvel_S, maxvel_L); if (maxvel_S > maxvel_L) fastest_p = fastest_p_S; else fastest_p = fastest_p_L; } if (topen>=t) // if topen is in the 'future' { // viscous damping during filling/settling stage: //if (ivdamp==1 && pe <= 0.5 initial_pe && ke > 1.0) if (ivdamp==1 && pe <= 0.5 initial_pe) { for (i=1; i<=np; i++) { cx[i] = vdc cx[i]; if (cy[i] > 1.0) cy[i] = vdc cy[i]; cz[i] = vdc cz[i]; wx[i] = vdc wx[i]; 202 PAGE 203 wy[i] = vdc wy[i]; wz[i] = vdc wz[i]; } } /*CRITERIA FOR OPENING HOPPER OUTLET */ if ( ke < 0.003 np + 0.5 && maxvel < 55.0 && t > 3000 && ifill!=2 ) // t>3000 to prevent opening with i.c.'s { iopen=1; writecount8 = 50000; // set writecount8 to a large number so that the initial state is recorded in filename_par_states.dat } /*PERFORM VARIOUS OPERATIONS AS HOPPER IS READY TO BE OPENED */ // Record level of free surface (c urrently only the yposition of th e uppermost particle is recorded) //if (iopen==1) { if (idiam==1) // record y position of uppermost particle { for (i=1; i<=np; i++) { if (y[i]>max_y) max_y=y[i]; } cout<<"Max Initial fill height is = "< PAGE 204 // Check mass fraction of initial samples on a spatial 2D grid ///* ARRAY DEFINITIONS for sample, p: sample_data[p][1] = x position of sample sample_data[p][2] = y position of sample sample_data[p][3] = z position of sample sample_data[p][4] = total number of particles in sample sample_data[p][5] = number of large particles in sample sample_data[p][6] = total mass of particles in sample sample_data[p][7] = mass of large particles in sample sample_data[p][8] = standard deviation of mass in sample sample_data[p][9] = mass fraction fines in sample, xi sample_data[p][10]= xi/xf */ for (i=0; i<=num_sam; i++) { for (j=0; j<=10; j++) { sample_data[i][j]=0.0; } } // take "samples" of initial fill, calculate mass fraction sample_data[1][1] = 0.5*(hdia2*sam_rad);// initial sample location sample_data[1][2] = sam_rad; sample_data[1][3] = 0.0; for (p=1; p<=num_sam; p++) { sam_loc=0; mass_ave = 0; if (p!=1) { sample_data[p][1] = sample_data[p1][1] + 0.25*(hdia2*sam_rad); sample_data[p][2] = sample_data[p1][2] ; sample_data[p][3] = 0.0; if (sample_data[p][1] >= 0.5 hdia) { sample_data[p][1] = 0.5*(hdia2*sam_rad); sample_data[p][2] += 0.25*(hdia2*sam_rad); } } for (i=1; i<=np; i++) { if ( sqrt( sq(x[i]sample_data[p][1 ]) + sq(y[i]sample_data[p][2]) ) < sam_rad ) { sample_data[p][4]++;//total count in sample sample_data[p][6] += mass[i];// total mass in sample 204 PAGE 205 if ( i <= nLG ) { sample_data[p][5]++;//count of large par ticles in sample should be OK for distributions particles in sample sample_data[p][7] += mass[i];// mass of large particles in sample } sam_loc++; if (sam_loc >= 9000) { cout <<"abort...sam_loc is too large for array sample_masses[9000]"< PAGE 206 { ofstream sample_out ("initial_samples.dat", ios::out); sample_out<<" x y z N N_LG M_tot M_LG S.D.mas*1000 xi xi/xf\n"; for (i=1; i<=num_sam; i++) { sample_out<<" "; for (j=1; j<=10; j++) { sample_out< PAGE 207 if (indy == 0) strcolormap < PAGE 208 min_pwoverlap_set = min_pwoverlap; max_pwoverlap_set = max_pwoverlap; sum_pwoverlap_set = sum_pwoverlap; part_collisions_set = part_collisions; total_wall_collisions_set = total_wall_collisions; overlap_warning = 0; min_overlap = 1; max_overlap = 0; sum_overlap = 0; min_pwoverlap = 1; max_pwoverlap = 0; sum_pwoverlap = 0; part_collisions = 0; total_wall_collisions = 0; pp_collisions_segn_stats = 0; pw_collisions_segn_stats = 0; // Finally, open hopper outlet //topen=t; // open slide gate // Set a counter for mass fraction output //mass_update_time = topen; writecount8=500000; // so temppos is written out when the slide gate first opens } } /*CRITERIA FOR ENDING RUN */ //OLD METHOD: a percentage of np //if ( (done == 0 && nrecyc!=0 && num_r ecycled > nrecyc*np ) (nrecyc==0 && num_discharge>=0.99*np && tfinal>t+5000 && hha <= 40)(nrecyc==0 && num_discharge>=0.92*np && tfinal>t+5000 && hha > 40 && hha <=70)  (nrecyc==0 && num_discharge>0.8*np && tfinal>t+5000 && hha > 70) ) // set endpoint of movie and simulation if (t PAGE 209 if ( (done == 0 && nrecyc!=0 && num_recycled > nrecy c*np )  (nrecyc==0 && tfinal>t+5000 && t prev_disch_time > 0.075/dt && t > topen) ) // set endpoint of movie and simulation { cout < PAGE 210 if (t>=topen && (topent)%10==0) // every 10 time steps after topen { for (i=1; i<=np; i++) // calculate average velocities { ave_velocity[i] += sqrt( sq(cx[i]) + sq(cy[i]) + sq(cz[i]) ); } num_velocities++; // increment counter if ( curr_randit==1 ) { for (i=1; i<=np; i++) { if (y[i] > 0) { xz_sq = sqrt( sq(x[i]) + sq(z[i]) ); k=m=0; // find k and m based on particle position if (ibc==3) // CAUTION: analysis in cylindrical bin not yet completed / working!!!!!! { do k++; while ( xz_sq > bin_sift[k][0][1]); // && xz_sq < bin_sift[k][0][1] ); } if (ibc==2) { do k++; while ( x[i] > bin_sift[k][0][1]); } do m++; while ( y[i] > bin_sift[0][m][1]); // && y[i] < bin_sift[0][m][1] ); /* cout << "i,x,y,z,cx,cy,cz = "<nLG) 210 PAGE 211 { bin_sift[k][m][7]++; bin_sift[k][m][8] += cx[i]; bin_sift[k][m][9] += cy[i]; // bin_sift[k][m][10] += cz[i]*sign(z[i]); bin_sift[k][m][11] += mass[i]; } } } } } writecount8++; if ( writecount8 >= twrite ) { writecount8 = 0; cout.setf(ios::showpoint); ofstream temp_out (filename12.c_str(), ios::out); for (i=1; i<=np; i++) { temp_out< PAGE 212 f_x = sqrt( sq(bin_sift[k][m][8]) + sq(bin_sift[k][m][10]) ) / bin_sift[k][m][7]; c_x = sqrt( sq(bin_sift[k][m][3]) + sq(bin_sift[k][m][5]) ) / bin_sift[k][m][2]; } f_y = bin_sift[k][m][9] / bin_sift[k][m][7]; c_y = bin_sift[k][m][4] / bin_sift[k][m][2]; sift_vel < PAGE 213 // clear bin_sift array for (k=1; k<=num_bin_x_sift+1; k++) { for (m=1; m<=num_bin_y_sift; m++) { for (w=2; w<=12; w++) bin_sift[k][m][w] = 0; } } } // close if curr_randit ==1 if (imovie==1) { if (t PAGE 214 partpos_out <<"0.0\t"<<0.2*hht<<"\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0"< PAGE 215 else ave_vel_out < PAGE 216 { if ( y[i] > 11.9 && y[i] < 12.1 ) vel_distn_out1 << t*dt<<" "< PAGE 218 flowdistn[distnlocation][5][curr_randit] = avedischvel / num_disch_in_massbin; flowdistn[distnlocation][6][curr_randit] = avedischmass / num_disch_in_massbin; flowdistn[distnlocation][7][curr_randit] = avedischdens / num_disch_in_massbin; flowdistn[distnlocation][8][curr_randit] = cummassdisch; flowdistn[distnlocation][10][curr_ randit] = pp_collisions_segn_stats; flowdistn[distnlocation][11][curr_randit] = pw_collisions_segn_stats; flowdistn[distnlocation][0][curr_randit] = 1.0; flowdistn[distnlocation][9][curr_randit] = avedischmassS / (avedischmassS + avedischmassL ) / x_fines; cout << "\nSegregation data...(for material left in bin) Iteration = "< PAGE 219 // if (idd == 2 && idiam == 1) flow_out< PAGE 220 time_rgb_out <<" 0 "<<(int)floor(0.5*255/fixed_time*time_k*5 + 255*2.5)<<" "<<(int)floor(255/fixed_time*time_k*5765)< PAGE 221 for (i=0; i<=13; i++) { sum_mass < PAGE 222 report_out << "\t\tMaximum pw settling phase over lap: \t" < PAGE 223 if (idiam==2) { report_out << "\n\t\t\t\t\t\tbut the maximum particle heights were:"< PAGE 224 report_out << "\t\tdia << dia < PAGE 225 { timedistn_out < PAGE 226 rtd_out << np << 0 100 100 100"< PAGE 227 } cout << "\nAbs(xVelocity) Distribution\n"; for (k=0; k<=num_bina; k++) { cout << setprecision(4) << binax[k][1] <<" "<< binax[k][2]< PAGE 228 ovrlp_out << setprecision(4) << bin_o[k][1] <<" "<< setprecision(2) << bin_o[k][2]< PAGE 229 cout << endl<< " "< PAGE 230 if (ibc==4) { cout << "\nThe walls underwent "< PAGE 231 cout << setw(7)< PAGE 232 cout << "ktw_f << ktw_f < PAGE 233 writecount6=count3=count8 = 0; for (i=1; i<=np; i++) skip_search[i] = skip_part[i] = 0; imovie=0; // turn off movie loops so excess files are not generated in subsequent runs } while (curr_randit < nreps+1); // do (entire program) while (more random iterations are called for) return(0); } void ParticleWallCollisions(double y_origin, double _hht, double _hdia, double _hod, double _hha, int open) { double _hhar, _tanhhar, _sinhhar, _coshhar; _hhar = _hha Pi / 180; _coshhar = cos(_hhar); _sinhhar = sin(_hhar); _tanhhar = tan(_hhar); if (ibc==2) // hopper { for (i=1; i<=np; i++) { n[3] = 0.0; if (skip_search[i]==0) { //Left bin wall if (x[i] <= 0.5 (_hdiadiam[i]) && y[i]>0+y_origin) { wall_num=8; overlap = 0.5*_hdia x[i] + 0.5*diam[i]; n[1] = 1.0; n[2] = 0.0; n_magn = 1.0; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } // Right Bin Wall if (x[i] >= 0.5 (_hdiadiam[i]) && y[i]>0+y_origin ) { wall_num=0; overlap = 0.5*diam[i](0.5*_hdiax[i]); n[1] = 1.0; n[2] = 0.0; n_magn = 1.0; 233 PAGE 234 ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } //Right hopper wall if ( ( _hha<90 && y[i]<1/_tanhhar 1.0*(x[i] 0.5*_hod) + 0.5*diam[i]/_sinhhar + y_origin && y[i] > 1.0*_tanhhar*(1.0*x[i] + 0.5*_hod) && y[i]>0+y_origin )  ( _hha==90.0 && y[i]<0.5*diam[i] && y[i]>y_origin && x[i]>0.5*_hod) ) { wall_num=1; overlap = 1.0*_sinhhar*(y[i](1/_tanhhar*(x[i]0.5*_hod)+0.5*diam[i]/_sinhhar)y_origin); n[1] = 1.0 _coshhar; n[2] = _sinhhar; n_magn = 1.0; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } // Left Hopper Wall if ( (_hha<90 && y[i] < 1.0/_t anhhar (1.0*x[i]0.5*_hod ) + 0.5*diam[i]/_sinhhar + y_origin && y[i] > 1.0*_tanhhar*(x[i]+0.5*_hod) && y[i]>0+y_origin)  (_hh a==90.0 && y[i]<0.5*diam[i] && y[i]>y_origin && x[i] < 0.5*_hod) ) { wall_num=7; overlap = 1.0*_sinh har*(y[i](1/_tanhhar*(x[i]0.5*_ho d)+0.5*diam[i]/_sinhhar)y_origin); n[1] = _coshhar; n[2] = _sinhhar; n_magn = 1.0; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } // slide gate (if hopper is NOT OPEN) if (open!=1) { if ( y[i] < 0.5*diam[i]+y_origin && x[i]<0.5*_hod && x[i]>0.5*_hod) { wall_num=4; overlap = 0.5*diam[i]y[i]+y_origin; n[1] = 0.0; n[2] = 1.0; n_magn = 1.0; ResetWallData(); 234 PAGE 235 NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } } // orifice walls and the points btw the orifice wa lls and the hopper walls (if the hopper is OPEN) if (open==1) { //point between right hopper wall and right orifice wall if ( ( _hha<90 && y[i]<1.0*_tanhhar*(x[i] 0.5*_hod) + y_origin && y[i]>0+y_origin && sq(0.5*_hodx[i])+sq(y[i]y_origin) < 0.25*diam[i]*diam[i] )  ( _hha==90.0 && x[i]<0.5*_hod && y[i]>0+y_origin && sq(0.5*_hodx[i])+sq(y[i]y_origin) < 0.25*diam[i]*diam[i] ) ) { wall_num = 2; overlap = 0.5*diam[i] sqrt( sq(0.5*_hodx[i]) + sq(y[i]y_origin)); n[1] = 0.5*_hod+x[i]; n[2] = y[i]y_origin; n_magn = sqrt ( sq(n[1]) + sq(n[2]) ); n[1] /= n_magn; n[2] /= n_magn; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } // The point between left hopper wall and left orifice wall if ( ( _hha<90 && y[i]<1.0*_tanhhar*(x[i] + 0.5*_hod) + y_origin && y[i]>0+y_origin && sq(0.5*_hodx[i])+sq(y[i]y_origin) < 0.25*diam[i]*diam[i] )  ( _hha==90.0 && x[i] > 0.5*_hod && sq(0.5*_hodx[i])+sq(y[i]y_origin) < 0.25*diam [i]*diam[i] && y[i]>0+y_origin ) ) { wall_num = 6; overlap = 0.5*diam[i] sqrt( sq(0.5*_hodx[i]) + sq(y[i]y_origin)); n[1] = 0.5*_hod+x[i]; n[2] = y[i]y_origin; n_magn = sqrt( sq(n[1]) + sq(n[2]) ); //sqrt( sq(0.5*_hod x[i]) + sq(y[i])); n[1] /= n_magn; n[2] /= n_magn; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } // right orifice wall if (y[i]<=0 && x[i]>=0.5*(_hoddiam[i]) && y[i]>=0.17*_hod0.025)// 'upper' hopper need not be included { wall_num=3; overlap = 0.5*_hod + x[i] + 0.5*diam[i]; 235 PAGE 236 n[1] = 1.0; n[2] = 0.0; n_magn = 1.0; // cout << "There has just been a collision btn particle "<=0.17*_hod0.025)// 'upper' hopper need not be included { wall_num = 5; overlap = 0.5*_hod x[i] + 0.5*diam[i]; n[1] = 1.0; n[2] = 0.0; n_magn = 1.0; // cout << "There has just been a collision btn particle "<= 0.0+y_origin && xz_sq >= 0.5 (_hdia diam[i]) ) { overlap = 0.5*diam[i] 0.5*_hdia + xz_sq; wall_num = 0; n[1] = 1.0 x[i]; n[2] = 0.0; n[3] = 1.0 z[i]; n_magn = sqrt( sq(n[1]) + sq(n[3]) ); // exclude n[2] always zero 236 PAGE 237 for (b=1; b<=3; b++) n[b] /= n_magn; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } // hopper walls if ( (_hha < 90 && y[i] < 1/_tanhhar*(xz_sq 0.5*_hod) + 0.5*diam[i]/_sinhhar+y_origin && y[i] > 1.0*xz_sq*_tanhhar + 0.5*_hod*_tanhhar && y[i]>0+y_origin)  (_hha==90 && y[i]<0.5*diam[i] && y[i]>0+y_origin && xz_sq>0.5*_hod) ) { //old (equivalent): overlap = 1.0*_sinhhar*(y[i] (1/_tanhhar*(xz_sq0.5*_hod) + 0.5*diam[i]/_sinhhar)y_origin); overlap = 0.5*diam[i] _sinhhar (0.5*_hod/_tanhhar (xz_sq/_tanhhar y[i]) ); wall_num = 1; n[1] = 1.0 x[i] _coshhar / xz_sq; n[2] = _sinhhar; n[3] = 1.0 z[i] _coshhar / xz_sq; //cout<<"\n t, i, x vector: "< PAGE 238 n_magn = sqrt( sq(n[1]) + sq(n[3]) ); // exclude n[2] always zero for (b=1; b<=3; b++) n[b] /= n_magn; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } //point between hopper wall and orifice wall if ( ( _hha<90 && y[i]<1.0*_tanhhar*(xz_sq 0.5*_hod) + y_origin && y[i]>0+y_origin && sq(0.5*_hodxz_sq)+sq(y[i]y_origin) < 0.25*diam[i]*diam[i] )  (_hha==90.0 && xz_sq<0.5*_hod && sq(0.5*_hodxz_sq)+sq(y[i]y_origin) < 0.25*diam[i]*diam[i] && y[i]>y_origin) ) { par = atan ( (y[i]y_origin) / (0.5*_hod xz_sq)); wall_num = 2; overlap = 0.5*diam[i] sqrt( sq( 0.5*_hod xz_sq ) + sq(y[i]y_origin) ); n[1] = 1.0 x[i] cos(par) / xz_sq; n[2] = sin(par); n[3] = 1.0 z[i] cos(par) / xz_sq; n_magn = sqrt ( sq(n[1]) + sq(n[2]) + sq(n[3]) ); // probably don't need this, n should already be normalized by the above def'ns for (b=1; b<=3; b++) n[b] /= n_magn; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } } // slide gate (if hopper is NOT OPEN) if ( open!=1 ) { if ( y[i] < 0.5*diam[i]+y_origin && xz_sq < 0.5*_hod) { overlap = 1.0 y[i] + 0.5*diam[i] + y_origin; wall_num = 4; n[1] = 0.0; n[2] = 1.0; n[3] = 0.0; n_magn = 1.0; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } } } // closes if pw_nn_list==1 } // close for i } // close if ibc==3 238 PAGE 239 if (ibc==4) // pharma hopper with ho pper section axis offset / tilted: { // ______ //   // ______ // \  // \  // \  // \__ // __ //phcm double _phir, _sin_phir, _cos_phir; double _psir; for (i=1; i<=np; i++) { if (skip_search[i]==0) { xz_sq = sqrt( sq(x[i]) + sq(z[i]) ); xz_varoff_sq = sqrt( sq( x[i] + _tanhhar*y[i]) + sq(z[i]) ); Hdist = sqrt( sq(z[i]) + sq( 1.0*alph/h2*y[i] x[i] ) ); gam = y[i]*_tanhhar + _hod; xz_off_sq = sqrt( sq(x[i]+alph) + sq(z[i]) ); // discharge chute if (y[i]>1.0*h1 && y[i]<0 && xz_sq>0.5*(_hoddiam[i])) { overlap = 0.5*diam[i] 0.5*_hod + xz_sq; wall_num = 3; // cout <<"\ncollision with discharge chute walls..."< PAGE 240 n[2] = 0.0; n[3] = 1.0 z[i]; n_magn = sqrt( sq(n[1]) + sq(n[3]) ); // exclude n[2] always 0.0 for (b=1; b<=3; b++) n[b] /= n_magn; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } // hopper walls if (y[i]>0.0 && y[i] 0.5*(gamdiam[i])) { nn = alph y[i] / h2; // distance from x=0,z=0 to hopper axis at this y[i] value _psir = atan(z[i]/(x[i]+nn)); // azimuthal angle (from RHS rotating forward to +z) if ( x[i]+nn<0 ) _psir += Pi; _phir = 0.5*_hhar ( cos(_psir) 1.0 ); // local hopper half angle (rad) _cos_phir = cos(_phir); _sin_phir = sin(_phir); overlap = 0.5*diam[i] 0.5*gam + Hdist; wall_num = 1; // cout <<"\ncollision with hopper wall..."<PAGE 241 241 n[1] = 0.0; n[2] = 1.0; n[3] = 0.0; n_magn = 1.0; ResetWallData(); NormalAndVelocity(i,0); SoftPWCollisionModel(i,wall_num); ParticleWallColl(i); } } } } } } PAGE 242 LIST OF REFERE NCES Ahn et al., 2008 H. Ahn, Z. Basaranoglu, M. Yilm az, A. Bugutekin and M.Z. Gul, Experimental investigation of gr anular flow through an orifice, Powder Technology, 186 (2008), pp. 6571. Anand et al., 2008 A. Anand, J.S. Curtis, C.R. Wassgren, B.C. Hancock and W.R. Ketterhagen, Predicting discharge dynamics fr om a rectangular hopper us ing the discrete element method (DEM), Chem. Eng. Sci. 63 (2008), pp. 58215830 Artega and Tuzun, 1990 P. Artega and U. Tuzun, Flow of binary mixtures of equaldensity granules in hopperssize segregation, flowing density and discharge rates., Chem. Eng. Sci., 45 (1990), pp. 205223. Barletta et al., 2003 D. Barletta, G. Donsi, G. Ferrari and M. Poletto, On the role and the origin of the gas pressure gradient in the discharge of fine solids from hoppers., Chem. Eng. Sci ., 5 (2003), pp. 52695278. Baxter et al., 2000 J. Baxter, H. AbouChakra, U. Tuzun, and B. M. Lamptey, A DEM simulation and experimental strategy for solving fine powder flow problems, Trans. Inst. Chem. Eng. 78(A) (2000), pp. 10191025. Beverloo et al., 1961 W.A. Beverloo, H.A. Leniger and J. va n de Velde, The flow of granular solids through orifices., Chem. Eng. Sci. 15 (1961), pp. 260269. Brown and Richards, 1960 R.L. Brown and J.C. Richards, Profile of flow of granules through apertures., Trans. Inst. Chem. Engrs. 38 (1960), pp. 243250. Brown and Richards, 1965 R.L. Brown and J.C.Richards, Kinematics of the flow of dry powders and bulk solids, Rheologica Acta, 4 (1965), p.153. Carson et al., 1986 J. W. Carson, T. A. Royal, and D. J. Goodwill, Understanding and eliminating particle segregation problems, Bulk Solids Handling 6 (1986), pp. 139144. 242 PAGE 243 Crewdson et al., 1977 B. J. Crewdson, A. L. Ormond and R. M. Nedderman, Airimpeded discharge of fine particles from a hopper., Powder Technology, 16, Issue 2 (1977), pp.197207. Cumberland and Crawford, 1987 D. Cumberland and R. Crawfo rd, The Packing of Particles, Handbook of Powder Technology., Vol. 6 (1987), Elsevier, New York Cundall and Strack, 1979 P.A. Cundall and O.D.L. Strack A discrete numerical model for granular assemblies., Gotechnique, 29 (1979), pp. 4765. Davies and Foye, 1991 C.E. Davies and J. Foye, Flow of granular material through vertical slots., Trans. Inst. Chem. Engrs., 69 (1991), pp. 369373. de Silva et al., 2000 S. de Silva, A. Dyry, and G. G. Enstad, Segregation mechanisms and their quantification using segregation testers. In IUTAM Symposium on Segregation in Granular Flows ed. by A. D. Rosato & D. L. Blackmore. Boston: Kluwer Academic Publishers (2000), 1129. Dias et al., 2004 R.P. Dias, J. A. Teixeira, M. G. Mota and A. I. Yelshin, Particulate binary mixtures: Dependence of packing porosity on particle size ratio., Ind. Eng. Chem. Res ., 43 (2004), pp. 79127919. Dziugys and Peters, 2001 A. Dziugys and B. Peters, A new approach to detect the contact of two dimensional elliptical particles, Int. J. Numer. Anal. Met ., 25(2001), pp. 14871500. Franklin and Johanson, 1955 F.C. Franklin and L.N. Johans on, Flow of granular material through a circular orifice., Chem. Eng. Sci., 4 (1955) (3), pp. 119129 Hermann and Luding, 1998 H. J. Herrmann and S. Luding, Modeling granular media on the computer, Continuum Mech. Thermodyn ., 10 (1998), pp. 189231. 243 PAGE 244 Hopkins and Tuhkuri, 1999 M. A. Hopkins and J. Tuhkuri, Compression of floating ice fields J. Geophys.Res., 104(C7) (1999), pp. 15,81515,825. Hsiau and Yang, 2003 S.S. Hsiau and S.C.Yang, Numerica l simulation of selfdiffusion and mixing in a vibrated granular bed with the cohesive effect of liquid bridges, Chem. Eng. Sci. 58 (2003), pp. 339351 Humby et al., 1998 S. Humby, U. Tzn and A.B. Yu, Prediction of hopper discharge rates of binary granular mixtures., Chem. Eng. Sci. 53 (1998), pp. 483494. Huntington and Rooney, 1971 A.P Huntington. N.M.and Rooney, Chemical Engineering Tripos part 2, Resear ch Project Report, University of Cambridge, 1971 Janssen, 1895 H.A Janssen, Versuche ber getreidedruck in silozellen, Z Ver.dt.Ing ., 39 (1895), p.1045. Johanson, 1996 J. R. Johanson, Predicting segregation of bimodal particle mixtures using the flow properties of bulk solids, Pharm. Technol ., 20 (1996), pp. 4657. Ketchum, 1929 M.S. Ketchum, The Design of walls, bi ns and grain elevators., McGrawHill, New York, 1929. Ketterhagen et al., 2005 W.R. Ketterhagen, J.S. Curtis a nd C.R. Wassgren, Stress results from twodimensional granular shear flow si mulations using various collision models., Phys. Rev., E. 71 (2005), p. 063307. Ketterhagen et al., 2008 W.R. Ketterhagen, J.S. Curtis C.R. Wassgren and B.C. Hancock, Modeling granular segregation in flow fr om quasithreedimensional, wedgeshaped hoppers, Powder Technology 179 (3) (2008), pp.126143. KruggelEmden et al., 2008 H. KruggelEmden, S. Wirtz, V. Scherer, A study on tangential 244 PAGE 245 force laws applicable to the discrete elem ent method (DEM) for materials with viscoelastic or plastic behavior, Chemical Engineering Science, 63 (2008) p. 1523. Lian et al., 1993 G. Lian, C. Thornton, M.J. Adams, A theoretical study of the liquid bridge forces between two rigid spherical bodies, J. Colloid Interface Sci 161 (1993), pp. 138147. Lian et al., 1998 G. Lian, C. Thornton, M.J. Adams, Discrete particle simulation of agglomerate impact coalescence, Chem. Eng. Sci. 53(19) (1998), pp. 33813391. Lin and Ng, 1997 X. Lin and T.T. Ng, A threedime nsional discrete element model using arrays of ellipsoids, Geotechnique 47(2) (1997), pp. 319. Luding et al., 2003 S. Luding, R. Tykhoniuk, and J. Toma s, Anisotropic material behavior in dense cohesivefrictional powders, Chem. Eng. Technol., 26(12) (2003), pp. 12291232 Mason and Clark, 1965 G. Mason and W.C. Clark, Liquid bridges between spheres. Chem Eng. Sci. 20 (1965), pp. 859. Matuttis et al., 2000 H. G. Matuttis, S. Luding, and H. J. Herrmann, Discrete element simulations of dense packings and heaps made of spherical and nons pherical particles, Powder Technology, 109 (2000), pp. 278. Matuttis and Schinner, 2001 H. G. Matuttis and A. Schinne r, Particle simulation of cohesive granular materials, Int. J. Mod. Phys. C, 12(7) (2001), pp. 10111021. McCarthy, 2003 J.J. McCarthy, Micromodeling of cohesive mixing processes, Powder Technology 138 (2003), pp. 63 67 Mehrotra and Sastry, 1980 V.P. Mehrotra, K.V.S. Sastry Pendular bond strength between unequal sized spherical particles, Powder Technology, 25 (1980), pp. 203214. Mikami et al., 1998 T. Mikami, H. Kamiya, and M. Horio, Numerical simulation of cohesive 245 PAGE 246 powder behavior in a fluidized bed, Chem. Eng. Sci ., 53 (1998), pp. 19271940 Myers and Sellers, 1971 M.E Myers and M. Sellers, Chemical Engineering, Tripos Part 2, Research Project Report, Un iversity of Cambridge, 1971. Nase et al., 2001 S.T. Nase, W.L. Vargas, A. A. Abatan and J. J. McCarthy, Discrete characterization tools for c ohesive granular material, Powder Technology 116 (2001), pp. 214223 Nedderman and Laohakul, 1980 R.M. Nedderman and C. Laohakul, The thickness of the shear zone of flowing granular materials., Powder Technology, 25 (1980), p. 91. Nedderman et al., 1982 R.M. Nedderman, U. Tzn, S.B. Savage and G.T. Houlsby, The flow of granular materials I di scharge rates from hoppers., Chem. Eng. Sci., 37 (1982) (11), pp. 15971609. Newton et al., 1945 R.H. Newton, G.S. Dunham G.S. and T.P. Simpson, Trans Am. Inst. Chem Engrs ., 41 (1945), p.215. Nguyen et al., 1979 T.V. Nguyen, C. Brennen and R.H. Sabersky, Gravity flow of granular materials in conical hoppers., Journal of Applied Mechanics, 46 (1979), pp. 529535. Ottino and Khakhar, 2000 J. M. Ottino and D. V. Khakhar, Mixing and segregation of granular materials, Annu. Rev. Fluid Mech ., 32 (2000), 5591. Popatov and Campbell, 1998 A. V. Potapov and C. S. Campbell, A fast model for the simulation of nonround particles, Granul. Matter 1(1998), pp. 914. Rausch, 1948 J.M.Rausch, Ph.D. Thesis, Princeton University, 1948. Ristow, 1997 G.H. Ristow, Outflow rate and wall stress for twodimensional hoppers., Physica A, 235 (1997), pp. 319326. 246 PAGE 247 Rose and Tanaka, 1959 H.E. Rose and T. Tanaka, Rate of discharge of granular materials from bins and hoppers., The Engineer 208 (1959), pp. 465469. Seville et al., 2000 J.P.K. Seville, C.D. Willett and P.C. Knight, Interparticle forces in fluidisation: a review, Powder Technology, 113 (2000), pp. 261268 Shi and McCarthy, 2008 D. Shi and J.J. McCarthy, Numerical simulation of liquid transfer between particles, Powder Technology 184 (2008), pp. 6475 Soulie et al., 2006 F. Soulie, F. Cherblanc, M. S. El Youssoufi and C. Saix, Influence of liquid bridges on the mechanical behaviour of polydisperse granular materials, Int. J. Numer. Anal. Meth. Geomech ., 30 (2006), pp. 213228 Spink and Nedderman, 1978 C. D. Spink and R. M. Nedderman, Gravity discharge rate of fine particles from hoppers., Powder Technology, 21, Issue 2 (1978), pp.245261. Suiker and Fleck, 2004 A.S.J. Suiker and N.A. Fleck Frictional Collapse of Granular Assemblies., J Appl. Mech ., 71 (2004), pp. 350358. Verghese and Nedderman, 1995 T. M. Verghese and R. M. Nedderman, The discharge of fine sands from conical hoppers., Chem. Eng. Sci., 50 (1995) (19), pp. 31433153. Walton and Braun, 1986 O.R. Walton and R.L. Braun, Vi scosity, granular temperature, and stress calculations for shearing assemblies of inelastic, frictional disks., J. Rheol., 30 (1986), pp. 949980. Walton and Braun, 1993 O. R. Walton and R. L. Braun, Simulation of rotarydrum and repose tests for frictional spheres and rigid s phere clusters, In Joint DOE/NSF Workshop on Flow of Particulates and Fluids (1993), pages 131148, Ithaca, NY. 247 PAGE 248 Weber et al., 2004 M. W. Weber, D. K. Hoffman, and C. M. Hrenya, Discreteparticle simulations of cohesive granular flow using a squarewell potential, Granul. Matter 6(4) (2004), pp. 239254. Willett et al., 2000 C.D. Willett, M.J. Adams, S.A. Johnson and J.P.K. Seville, Capillary bridges between two spherical bodies, Langmuir, 16 (2000), pp. 93969405. Williams and Pentland, 1992 J. R. Williams and A. P. Pentland, Superquadrics and model dynamics for discrete elements in interactive design, Eng. Comm., 9 (1992), pp. 115127. Yang et al., 2003 R. Y. Yang, R. P. Zou, and A. B. Yu, Numerical study of the packing of wet coarse uniform Spheres, AIChE Journal, 49(7) (2003), pp. 16561666. Zhou et al., 2005 Z. Zhou, H. Zhu, A. Yu, B. Wright, D. Pinson and P. Zulli, Discrete particle simulation of solid flow in a model blast furnace., ISIJ Int. 45 (2005), pp. 18281837. 248 PAGE 249 249 BIOGRAPHICAL SKETCH Anshu Anand was born in 1982 in Ranchi, India, the son of parents Surendra and Pushpa Singh. He attended the Indian Institute of Tec hnology, Bombay and earned a bachelors degree in Chemical Engineering in May 2005. Research experience gained during this time included internship at the Max Planck Institut fur Chemie Mainz Germany. The author then enrolled at University of Florida, Gainesville to pursue graduate studies in Chemical Engineering. While here, his industrial experience was expanded by an internship with Pfizer, during which time cohesive forces between particles were meas ured using the Atomic Force Microscopy. He received his PhD from the University of Florida in the summer of 2009. 