' ("WE—W" ‘,_7~.'- — v c--.» 3 A N... .4 . ' “:7 . I? H b ‘ I 4 :git‘iéii:* ~:¥X ‘1; 9.33;: 5 “~xs L..-‘ ,‘, A - 9 x‘tré a- -7, ‘ K\ #7} 1"0.' I ’4' ,- .J .5". This is to certify that the dissertation entitled FINITE ELEMENT SIMULATIONS OF CERAMIC POWDER COMPACTION AND SINTERING IN THE MAKING OF A MICRO HEAT EXCHANGER presented by CHEE KUANG KOK has been accepted towards fulfillment of the requirements for the Ph. D. degree in ENGINEERING MECHANICS IN MECHANICAL ENGINEERING 24%? L, 44/ T Major Professor’s Signature‘ flame Date MSU is an Affirmative Action/Equal Opportunity Institution "—-O-I-O-O-0-.-.-0-0-O-O-I-0-0-O-O—O-O-O-U-I-.-.-.-.-O-0-0-U-O-C-O-I-O-O-O-C-O-C-.—.-.-I-I’. - . ,-..._ LIBRARY Michigan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE II"' ‘0 ‘ 3: 6/01 cJCIFIC/DateDuepGS—sz FINITE ELEMENT SIMULATIONS OF CERAMIC POWDER COMPACTION AND SINTERING IN THE MAKING OF A MICRO HEAT EXCHANGER By CHEE KUANG KOK A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPH'-' Engineering Mechanics in Mechanical Engineering 2004 ABSTRACT FINITE ELEMENT SIMULATIONS OF CERAMIC POWDER COMPACTION AND SINTERING IN THE MAKING OF A MICRO HEAT EXCHANGER By CHEE KUANG KOK , I kOne major problem in processing ceramics using conventional powder- based compaction and sintering technique is the non-uniform relative density distribution and the residual stresses induced during these processes. The non- uniforrnities may cause high differential sintering rate, introducing cracks and shape distortion in the sintered product. The capability in predicting the relative density distribution and residual stresses will therefore enhance our ability to design and fabricate ceramic components with internal cavities?) In this thesis, finite element models for compaction process are developed based on theory of plasticity using the modified Drucker-Prager/cap model in soil mechanics and a recent hyperbolic cap model. A finite element model for sintering process based on viscoplasticity constitutive law is also presented. The models allow simulations to be performed before experiments to select feasible processes. The results from finite element simulation are consistent with findings of other researchers and in agreement with our experimental observations in making various features in our meso-scale heat-exchanger. The compaction models indicate that severe local shear yielding during unloading may have caused micro-cracks since the sites that yielded in shear during top punch removal are coincident with the sites where micro-cracks originate. The sintering model implies that sintered products are more uniform in their final relative density and residual stress distribution regardless of their geometry and their base states as inherited from the two different compaction models. In addition, simulations reveal that local effects such as differences in relative density, residual stresses and geometry do not seem to significantly affect the final grain size and sintering time needed to density a ceramic compact. An in-depth study of shear yielding and its relationship with compact cracking is needed to improve the current compaction model. To fine-tune the sintering model, sinter forging needs to be conducted to obtain experimental material parameters. Copyright by CHEE-KUANG KOK 2004 To my niece, Esther ACKNOWLEDGMENTS This material is based on research sponsored by the Air Force Office Research Laboratory, under agreement F49620-02-1-0025. The U. S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. I would like to express my gratitude for the financial support. In addition, I would like to thank my advisor, Dr. Patrick Kwon, for his trust and his vast support in this work. I could not have gone as far without his true understanding of my need for idea incubation time in making each baby step of progress in this work. It has been a tremendous pleasure to work for him. Dr. Eldon Case, Dr. HungYu Tsai, Dr. Ronald Averill, Dr. Farhang Pourboughrat, and my colleagues, YuHui Wang and YuWei Chi, have at various stages of this work made their timely contributions, and their efforts are greatly appreciated. Last but not least, I owe this work to all who have given me lasting hope, particularly my family and my Lord Jesus Christ. vi TABLE OF CONTENTS LIST OF TABLES ................................................................................... x LIST OF FIGURES ................................................................................ xi CHAPTER1 INTRODUCTION .................................................................. 1 1.1 Introduction ............................................................................. 1 .1. 1 Powder Compaction .................................................... 2 .1.2 Compact Sintering ...................................................... 5 /1./2 Background and Motivation ........................................................ 6 /1/ 3 Literature Review of Compaction Modeling ................................... 8 1 ..3 1 Constitutive Models in Literature for Compaction .............. 9 1.3.1.1 First-Order Model ........................................ 11 1.3.1.2 Elasticity Model ........................................... 12 1.3.1.3 Plasticity Model ........................................... 13 __ 1.3.2 Numerical Methods ................................................... 19 v1.4 Literature Review of Sintering Model .......................................... 21 1.4.1 Sintering Physics ...................................................... 21 1.4.2 Sintering Kinetics ...................................................... 22 1.4.3 Sintering Models ....................................................... 25 1.4.3.1 Background....-. ........................................... 25 1.4.3.2 Porous Viscosity Model .................................. 27 1.4.3.3 Sintering Potential ........................................ 29 1.4.4 Grain Growth Models ................................................. 33 1.5 Thesis Overview .................................................................... 35 CHAPTER 2 CONSTITUTIVE MODEL AND NUMERICAL SCHEME FOR COMPACTION SIMULATION .............................................. 36 2.1 Constitutive Model .................................................................. 36 2.1.1 Motivation ................................................................ 37 2.1.2 Hyperbolic Cap Model ................................................ 40 2.1.2.1 Linear Drucker-Prager ................................. 41 2.1.2.2 Cap Surface .............................................. 44 2.1.3 Modified Drucker-Prager/Cap Model ............................. 47 2.1.4 Differences Between Hyperbolic Cap Model and Modified Drucker-Prager/Cap Model ........................... 49 2.1.5 Plasticity ................................................................. 50 2.2 Numerical Scheme ................................................................ 53 2.2.1 Implementation ........................................................ 53 2.2.2 Yielding Algorithm .................................................... 55 2.2.2.1 Yielding on Linear Drucker-Prager Shear Surface ............................................ 55 vii 2.2.2.2 Yielding on Modified Hyperbolic Cap ............... 56 2.2.3 Linearization Moduli ................................................ 60 2.2.3.1 Linearization Moduli for Shear Yield Surface......61 2.2.3.2 Linearization Moduli for Hyperbolic Cap Surface ...................................................... 62 2.3 Implementation ..................................................................... 64 2.3.1 Adopted Zirconia Data ............................................... 64 .,2.3.2 Friction Between Powder Bed and Die Wall .................... 66 {2.3.3 Finite Element Approach ............................................. 67 CHAPTER 3 COMPACTION SIMULATION RESULTS ................................. 70 3.1 Verification ............................................................................. 70 3.1.1 Results and Discussions ............................................. 70 3.2 Compaction Without Graphite Piece .......................................... 72 3.2.1 Results and Discussions ............................................ 72 3.3 Compaction With Graphite Piece .............................................. 75 3.3.1 Results and Discussions ............................................ 77 3.3.1.1 Three Dimensional Compaction ...................... 77 3.3.1.2 Two Dimensional Compaction ........................ 79 3.3.1.2.1 Size and Orientation Effects ............... 79 3.3.1.2.2 Load Effects .................................... 83 3.4 Experimental Observations ...................................................... 84 3.4.1 Size and Orientation Effects ........................................ 85 3.4.2 Load Effects ............................................................. 86 3.5 Micro-crack Site Predictions ..................................................... 87 3.5.1 Unloading UniaxialIy-Pressed Sample ........................... 90 3.5.2 Shear Yielding and Crack Sites .................................... 91 3.6 Conclusions ......................................................................... 94 CHAPTER 4 CONSTITUTIVE MODEL AND NUMERICAL SCHEME FOR SINTERING SIMULATION ................................................... 96 4.1 Constitutive Model ................................................................. 96 4.1.1 Constitutive Model Background ................................... 96 4.1.2 Newtonian Fluid Model Basics .................................... 99 4.1.3 Newtonian Fluid Model to Sintering Model ................... 102 4.1.4 Shear and Bulk Viscosities ....................................... 104 4.1.5 Grain Growth Model ................................................ 110 4.2 Finite Element Model ............................................................ 111 CHAPTER 5 SINTERING SIMULATION RESULTS .................................. 113 5.1 Curve-Fit and Verification ....................................................... 113 5.2 Sintering of Powder Compact ................................................. 118 5.2.1 Coefficient of Effective Diffusion ....................... 118 5.3 Results and Discussions ....................................................... 119 5.3.1 Shrinkage and Shape Distortion ...................... 120 5.3.2 Sintering Time ............................................. 126 viii 5.4 Conclusions ........................................................................ 127 CHAPTER 6 CONCLUSIONS ............................................................... 128 6.1 Conclusions .......................................................................... 128 , 6.2 Future Works ....................................................................... 129 APPENDIX A: Partial Derivatives for Hyperbolic Cap Surface ...................... 131 APPENDIX B: Newton Iteration for Compaction Constitutive Model ............... 132 APPENDIX C1: Linearization Moduli for Shear Yield Surface ....................... 135 APPENDIX C2: Linearization Moduli for Hyperbolic Cap Surface .................. 138 APPENDIX D: Constitutive Laws for lncompressible Fluid ........................... 140 BIBLIOGRAPHY ................................................................................. 141 ix LIST OF TABLES Table 2.1 Mechanical properties of bulk zirconia ......................................... 64 Table 2.2 Geometry of simulated compacts ............................................... 68 Table 4.1 Particle sizes of 3 mol percent yittria-stabilized zirconia .................. 97 Table 4.2 Zirconia sintering data ............................................................. 108 Table 5.1 Average relative density and grain size for different cases ............. 125 LIST OF FIGURES Figure 2.1 Yield locus of the linear Drucker—Prager function .......................... 41 Figure 2.2 Hyperbolic-cap and linear Drucker-Prager yield functions .............. 44 Figure 2.3 Schematic of flow for the modified Drucker-Prager/cap model in p-q space ................................................................ 47 Figure 2.4 Schematic of flow for the linear Drucker/Prager model in p-q space ......................................................................... 55 Figure 3.1 Relative density distribution of the ejected right-half axisymmetric disc of height = 15.746mm and radius = 6.55mm from (A) Experimental results by Kim et al [9], (B): the hyperbolic model; (C): the modified Drucker-Prager/cap model. The disc is subjected to 100 MPa axial loading ........................................... 71 Figure 3.2 Relative density distribution of the ejected right-half axisymmetric disc of height = 15.746 mm and radius = 6.55 mm from (A): the modified Drucker-Prager/cap model and (B): the hyperbolic model. The disc is subjected to 23 MPa axial loading ........................................................................ 73 Figure 3.3 Relative density distribution of the ejected right-half axisymmetric disc of height: 4.114 mm and radius = 11.125 mm from (A): the modified Drucker-Prager/cap model and (B): the hyperbolic model. The cylinder is subjected to 23 MPa axial loading. Relative density value markers are the same as those in Figure 3.2 ............................................................. 74 Figure 3.4 Residual Mises stress distribution of the ejected right-half axisymmetric disc in Figure 3.3 (A): the DP model and (B): the hyperbolic model ........................................................ 74 Figure 3.5 Relative density distribution after ejection for a quarterly powder bed with embedded graphite piece subjected to 11.5 MPa. The empty spot is where the graphite resides. Powder bed initial thickness is 6.0 mm, and its radius 11.125 mm. The graphite piece has a cross-section of 1mm X 1m ........................ 78 Figure 3.6 Relative density distribution after ejection simulated using the modified Drucker-Prager/cap model (top row) and the hyperbolic model (bottom row). A graphite piece is xi embedded in different orientations and locations as in (A), (B) and (C). The graphite piece has a cross-section of 2 mm X 2 mm. All compacts are subjected to 23 MPa. The preform has a thickness of 6 mm .............................................. 80 Figure 3.7 Relative density distribution after ejection simulated using the modified Drucker-Prager/cap model (top row) and the hyperbolic model (bottom row). A graphite piece is embedded in different orientations and locations as in (A), (B) and (C). The graphite piece has a cross-section of 1 mm X 1 mm. All compacts are subjected to 23 MPa. The preform has a thickness of 6 mm ............................................. 81 Figure 3.8 Relative density distribution after ejection simulated using the modified Drucker-Prager/cap model (top row) and the hyperbolic model (bottom row). A graphite piece is embedded in different orientations and locations as in (A), (B) and (C). The graphite piece has a cross-section of 1 mm X 1 mm. All compacts are subjected to 11.5 MPa. The preform has a thickness of 6 mm ............................................. 84 Figure 3.9 Micro-cracks observed in both samples with graphite piece of cross-section 2 mm X 2 mm subjected to 23 MPa compaction load. Graphite piece is oriented flat in (A) and diagonally in (B) .................................................................................. 85 Figure 3.10 Micro-cracks observed in both samples with graphite piece of cross-section 2 mm X 2 mm subjected to 23 MPa compaction load. Graphite piece is oriented flat in (A) and diagonally in (B) ................................................................................ 86 Figure 3.11 Lowering the compaction pressure helps reducing or even eliminating micro-cracks in samples ......................................... 87 Figure 3.12 Loading and unloading paths the unaxiaIIy-pressed sample, 2 being the axial direction ...................................................... 90 Figure 3.13 Simulation results identifying regions that have yielded by shear (lighter area) before the ejection step under 11.5 MPa pressing with embedded graphite piece of the two sizes mentioned above ................................................... 92 Figure 3.14 Simulation results identifying regions that have yielded by shear (lighter area) before the ejection step under 23 MPa pressing with embedded graphite piece of the two sizes mentioned above ................................................... 93 xii Figure 5.1 Applying equal temperature to all nodes in the element ................ 114 Figure 5.2 Densification observed from experiment and from simulation of a unit cell using coefficients 02:02 and 04:03. ..................... 114 Figure 5.3 Scanning electron images depicting increasing grain size with respect to sintering time at 1475 degree Celsius ......................... 116 Figure 5.4 Grain growth observed from experiment using linear intercept method and from simulation of a unit cell .................................. 117 Figure 5.5 Grain growth observed from experiment using linear intercept method and from simulation of a unit cell, with sintering time plotted on log scale ............................................................... 118 Figure 5.6 Deformed configuration of sintered disc (solid lines) superimposed on the green compact (dotted lines) ..................... 120 Figure 5.7 Deformed configurations of sintered compacts with channels (solid lines) superimposed on the green compacts (dotted lines) ...................................................................... 121 Figure 5.8 Relative density distribution of compacted and finally sintered disc (not drawn to scale) ....................................................... 122 Figure 5.9 Hydrostatic stress distribution of compacted and finally sintered disc (not drawn to scale) ....................................................... 123 Figure 5.10 Relative density distribution of compacted and finally sintered disc ................................................................................. 124 xiii Chapter 1 INTRODUCTION 1.1 INTRODUCTION QAodem engineering ceramics such as alumina, zirconia, silicon carbide, tungsten carbide and etc. are making more usages in aerospace, automotive and electronics industry. Compared to metals, most covalent-bonded ceramics are harder, stronger, and possess high melting temperature, relatively low thermal conductivity, low fracture toughness and low breaking energy. They also exhibit low symmetry of crystal structure and limited or inadequate slip systems for plastic deformation. However, these properties of ceramic materials also make them more difficult to process. In this research, a meso-scale heat exchanger using ceramic powder is to be fabricated. This device contains a complex internal network of channels to circulate cooling fluids. The processing routine has been developed, where ceramic powder is poured into a mold and subjected to uniaxial compaction. lntemal meso-scale channels are introduced in the bulk of ceramic materials by embedding a fugitive phase on the surface of powder bed or inside the ceramic powder during compaction. When the compact is sintered to yield the final product, the fugitive phase burns away during sintering, leaving behind the channels in the shape of fugitive phase in the bulk of ceramics. Alternatively, the materials with the surface channels can be joined to form internal channels. The major drawback in this fabrication technique is the limitation in the channel shapes it can produce. As the shape of the channels becomes complex, powder flow around the complex fugitive phase becomes increasingly limited. This may lead to large crack formation around the fugitive phase in the final ceramic parts. The main goal of this thesis is to develop a computer model that simulates the compaction and sintering processes. The effects of powder compaction and powder flow may be studied from simulation results of density and residual stress distribution. 1.1.1 POWDER COMPACTION . Dry cold compaction is a common process for ceramic components that involves relatively simple technique to compact ceramic powder before sintering. However, the friction on the die wall and the design of mold often causes density variation and residual stresses within the compact. Heterogeneity in density results in non-uniform shrinkage of the powder compact during sintering, causing unpredictable dimensions of the final product. Coupled with high residual stress, regions of high density-gradient may lead to defects such as micro-cracks and capping effects. Tight tolerance on the final product dimension is essential for ceramic products of high quality. Defects in the final products are unacceptable and must be eliminated. Thus the ability to predict the density variation and residual stresses induced during the compaction process is highly desirable.) The compaction process can be described using numerical solution coupled with continuum mechanics based on the developments made in plasticity and soil mechanics for determination of stress distribution under complex loadings. The earliest plasticity models include the Extended Tresca criterion, the Extended Von Mises criterion (also known as the linear Drucker-Prager model) [1] and the Mohr—Coulomb criterion [2]. However, these models mostly capture soil behavior under shear loadings when soil is assumed to exhibit plastic behavior. The hardening of soil under hydrostatic stress has not been accounted for in these models. To capture the hardening effect, additional yield criterion has been introduced into later models. The Critical State model [3] and the Cam-Clay model [4] were developed. These models generally incorporate yield criteria for the shear loadings and cap yield criteria for triaxial loadings. The modified Drucker-Prager/cap model [5] is one of the cap models, which has been widely used in modeling soil compaction, and is available in the function library of commercial finite element code, such as ABAQUS. Researches in powder metallurgy have tried applying principles derived from soil particles to understand the behavior of metal and ceramic particles. Despite the similarities, there are some key differences among metal, soil and ceramic particles, as summarized by Bortzmeyer [6];. 1) 2) 3) Metal particles are typically large (about 100 um); and their packing lies between ‘loose’ and dense’. Since they are usually subjected to compaction deformation, their major microscopic compaction mechanism is the plastic deformation of the particles. Soil particles are often large (sand), arranged in loose or dense packing. They are usually subjected to shear deformation, and since the particles are rigid, the microscopic mechanism is the rearrangement of the particles. Ceramic particles are small (often less than 1 pm) and rigid. The van der Waals attraction among them is far greater than their weight, so that their relative density as poured into a mold may be as small as 20% (far less than loose packing). Like the soil particles, the microscopic mechanism is the rearrangement of particles. Bortzmeyer [6] suggested that ceramic powders should be treated as a distinct class of material that requires specific models or major modification of classical soil models. Even so, due to similarity in the mechanical behavior of ceramic powder and soil, modeling approaches by many recent researchers such as Aydin et al. [7,8], Kim et al. [9] and Zipse [10] are heavily based on soil mechanics models. 1 .1 .2 COMPACT SINTERING In modeling the sintering process, many researchers begin by describing its physical mechanism. The densification kinetics of a powder aggregate was first derived from the sintering of two contacting spheres by viscous flow or diffusion phenomenon. Analytical equations are proposed to investigate the two-sphere problem and the macroscopic behavior then derived. Because of the many variables affecting sintering, analytical solutions are only possible by making considerable geometric and diffusion flow field approximations rarely realized in practice. Barsown [11] observed that all earlier sintering models have these common traits: a) An assumed representative particle shape b) Surface curvature calculated as a function of geometry 0) Flux equation that depends on the rate-limiting (the slowest species diffusing along its fastest path) process. (I) Integrated of flux equation to predict the rate of geometric change. Numerical methods that give a more realistic description of the inter-particle neck geometry, the coupling of several diffusion mechanisms and the processing of a large randomly packed particles are developed thereon. These micro-structural models usually require empirical adjustments of the physical parameters for different materials. A lot of groundwork has been laid out in studying the relationship between the packing geometry and the relative density of metallic powder, and how diffusion mechanisms play a role in affecting the formers. These efforts include Artz [12] and Fischmeister et al. [13, 14], who first described the arrangement of particle centers by a radial distribution function, assuming spherical monosize powder particles. In their study, the powder compact densifies in an assumed particle radial growth mechanism, which leads to overlapping regions of neighboring particles, just like particles forming necks in physical observation. This overlapping volume is then distributed to the void space, assuming the dominant redistribution mechanism in operation [15]. These thought experiments have allowed for the quantitative description of densification rates according to some random packing geometry while taking account of the sintering mechanisms. Following their investigation, sintering diagrams [12] are later developed to highlight the dominant sintering mechanism given an arbitrary combination of sintering conditions. These efforts are further expanded by Li and Funkenbush [16,17] and Nair et al. [18, 19] to account for bimodal or unequal- sized hard powder particles. 1.2 BACKGROUND AND MOTIVATION A cost-efficient and accurate approach to manufacture a meso—scale ceramic heat exchanger is pursued in the thesis. Computer predictions of compaction and sintering in the fabrication of the micro heat exchanger may be important in designing the manufacturing processes and the product itself. Simulation results often provide insight into the ceramic material behavior under compaction and sintering. In addition to predicting density gradients, the numerical simulation provides insight into the evolution of density within the compact during the entire compaction process, which is difficult, if not impossible, to accomplish experimentally. The predicted evolution of the relative density and stress distributions within a compact during compaction can provide information on the likelihood of defects or fracture during ejection and subsequent sintering. In short, an accurate and efficient simulation model saves both the time and cost related to numerous testing of physical processes. In developing and implementing the computer model, this thesis is concerned mostly with the compaction and ejection processes, the subsequent sintering and the net-shape predictions of the final product. Other important elements such as powder flow for die filling and pre-sintered state machining will be considered in future research. The following Section 1.3 and Section 1.4 present literature reviews of the modeling works by previous researchers in compaction and sintering respectively. 1.3 LITERATURE REVIEW OF COMPACTION MODELING iCompaction modeling provides a basis for machine-loading specifications, and guidelines to help minimize defects such as warpage and visible cracks that form at the top of the greens compact, and micro-cracks that tend to occur at regions of high residual stress. The value and the predictive capability of any model are entirely dependent on the availability and the reliability of the required material input data, and the insight for how that material behaves under the considered loading conditions. Model validation is important, since ceramic powder compact behavior can be very complex to model accurately without relying on experimental input. This usually requires experimentally characterizing the compact and the compaction process. In choosing suitable constitutive models, the responses of the material to the environment and the applied stress states must be considered. They include: 1) The yielding, hardening, and failure characteristics of the frictional and compressible powder during compaction 2) The large strain formulation of the powder bed under compaction 3) The friction at the interface of the powder and the wall. 4) The elastic and possibly the plastic recovery within the compact upon completion of the compaction process. An analysis of the powder compaction process is therefore critically dependent on the validity of the chosen yield criteria and also the stress-strain relationships used to describe the behavior of the powder material after yielding. Aydin et al. [20], in a review on compaction models in literature, classifies the models into the following hierarchies: 1) First order models essentially based on continuum assumptions. 2) FEM models that are two or three-dimensional and are normally continuum based. 3) Distinct or discrete models that consider particle-particle and particle-wall interactions. 4) Boundary element models that appropriately handle large strains generated in certain compaction processes. All models require material input data for both the wall, the bulk powder, and the compact rheology, with more refined models requiring more input data. 1.3.1. CONSTITUTIVE MODELS IN LITERATURE FOR COMPACTION .«*‘”Granular material such as powder is known to show a great variety of behaviors L, when subjected to different conditions. Following the argument of Chen [21], real- life powder pouring and subsequent compaction processes may introduce some weak seams or loose spots in a ceramic compact. This is difficult to detect and its implication cannot be easily captured in simple theories. Even the most sophisticated theory may not describe fully the details of an actual process due to the imperfectness of the laboratory investigations. Simplifications and idealizations are therefore necessary to produce simple but accurate models that can represent the key features of the material at hand. These models, therefore, should not be expected to be valid over a wide range of conditions. Application of different models should be consistent with their assumptions. Following Chen’s proposal, the following three basic criteria for model evaluation may be considered: 1) Theoretical evaluation of the models with respect to the basic continuum mechanics principles to ensure consistency with the theoretical requirement of continuity, stability and uniqueness. 2) Experimental evaluation of the models with respect to their suitability in fitting experimental data from a variety of available tests, and the ease of determination of the material parameters from standard test data. 3) Numerical and computational evaluation of the models with respect to their implementation by computer algorithms and available hardware and software. 10 In a nutshell, when evaluating a model, a balance should be considered among the requirements in the continuum mechanics, the requirements of realistic representation of soil behavior during the experimental testing, and the requirements for the simplicity in application from a computational viewpoint. 1.3.1 .1 FIRST-ORDER MODEL The classical first-order model in compaction analysis, or the Jansen-Walker model [22], relates the stress transmission ratio of a cylindrical powder compact uniaxially pressed to the aspect ratio of the compact. The influence of powder-die wall friction and the radial and axial stresses are captured in this model. It is further modified by Briscoe and Ozkan [23] to incorporate a relationship between the powder-compact density and the compaction pressure. Since the work of Jansen-Walker, there have been more than 20 different modifications proposed to describe uniaxial compaction [20]. However, it seems that there is not a universal equation that fits into all loading conditions, since each model supposes different operative compaction mechanisms over the distinct range of applied stresses. An analytical equation based on fundamental understanding of powder compaction, which encompasses the knowledge of intrinsic flow properties of powder under stress, the local strength and stress distribution within a compact and the traction at particle-to—particle contact points is yet to be developed. 11 1.3.1.2 ELASTICITY MODEL A first step to break away from the limitation of first-order model is the use of elasticity model. Under a relatively short loading period, ceramic powder compact behavior may be idealized as time-independent. Elastic models are generally classified into two classes: 1) hyperelastic and 2) hypoelastic. The minimum requirement for a material to qualify as elastic is the existence of a one-to—one correspondence between stress and strain increments. The body returns to its original state of deformation when all the stresses are relaxed. Elastic materials satisfying this reversibility in the infinitesimal sense is said to be hypoelastic. More rigorously however, an elastic material that also satisfies the energy equation of thermodynamics is hyperelastic [24]. Aydin et al. [25] adopted a hypo-elastic nonlinear model in an attempt to simulate ceramic compact behavior. The material properties that define the stiffness matrix, Q" as in dg=_C_“’dg"’, include Young’s Modulus of Elasticity and Poisson’s ratio as a function of strain invariants. This model may be considered appropriate if the primary objective of the analysis is to predict the final state of the compact under compaction pressure but not the microscopic origins of the flow behavior of the powder. The simulation results do not compare well with the experimental results quantitatively, though the general trend of relative density 12 distribution is well captured. Aydin et al. [25] acknowledge the inability of the model to predict characteristic high and low-density regions that are observed experimentally in the lower and upper axial regions respectively in a cylindrical compact. It is believed that the disparity arises from the internal frictional effects, namely elements of radial flow must occur during compression, which is in part described by the internal angle of friction of a powder and is not modeled through the hypo-elastic approach. In soil mechanics, the linear theory of elasticity is the most commonly used constitutive model for soils, despite its shortcomings. According to Chen [24], the linear elastic model can be significantly improved by assuming a bilinear or higher order polynomial type of non-linear fit for the stress-strain relationship of soil in the form of secant formulations. This is known as the Cauchy type of elastic formulation. This type of model must be combined with criteria defining ‘failure’ of the soil. However, these formulations fail to identify plastic deformation when unloading occurs. Introducing loading criteria as in the deformation theory of plasticity can rectify this inadequacy to some extent. 1.3.1.3 PLASTICITY MODEL Adopting similar principles from soil mechanics, ceramic compact behavior has been commonly idealized as elastic-plastic. Owing to the differences between the 13 properties of ceramic powder as granular material and those of the idealized bodies, the classical plasticity theories must be modified and extended so that the unique preperties of ceramic powder compact are taken into consideration. Chen [24] claimed that the flow theory of plasticity was necessary to correctly account for the permanent plastic strain in addition to the elastic strain to model granular material behavior. The stress-strain law for a plastic material is essentially a relation between the stress and plastic strain increments in the current stress/strain state. This incremental relation is generally assumed to be homogenous and linear, thus a time-independent idealization that precludes viscous effects. The development of the incremental stress-strain relations in the flow theory of plasticity is based on three fundamental assumptions [26]: 1) The existence of initial and subsequent yield (loading) surfaces, f = 0 2) The formulation of a hardening rule that describes the evolution of subsequent loading surfaces, 3) A flow rule that specifies the general form of the plastic stress-strain relation. For plasticity models, an orthodox approach for formulating the flow rule is based on the assumption of separating the net strain increment into its elastic and plastic components. The plastic strain increment is then assumed to be 14 orthogonal to the stress gradient of the plastic potential function g as in daP’ = (123—. 80' For some materials, the plastic potential function is selected to be the yield function. In this case, the flow rule is associative and it satisfies the normality condition. However, for a granular material, this condition may lead to excessive volume change during shear deformation [2]. Thus a non-associative flow rule is usually required to model powder material shear behavior. Shima and Mimura [27] conducted rigorous testing on the silicon nitride ceramic powder and found that the normality of the strain rate vector to the yield surface they proposed almost holds, except for strain ratio exceeding certain limit. Bortzmeyer [6], however, claimed that the flow rule for granular materials is not generally associated. Non-associativity, he believed, is a direct consequence of the fact that the particles are rigid. This prevents overlapping of particles and as a result, Hill’s principle of maximum work needs to be modified to account for the constrained motion. Bortzmeyer [6] noted that the flow rule is slightly non- associated along the cap, and non-associativity becomes more severe as the stress state becomes more deviatoric. Along the critical state line, the law has to be non-associative if the material volume is constant during excessive shear deformation. 15 Bortzmeyer [6] also observed some irreversible deformation during unloading or reloading that produces the well-known hysteresis loops in the stress-strain relationship. This is not very important in ceramic technology, since the production requirements usually imply monotonic loadings. Caution has to be taken in neglecting the phenomenon because sample rupture during unloading or ejection may probably be related to this. Generally, two types of plasticity models are available in the literature, namely the single-yield function models and the multi-yield function models. Single-yield function model includes the well-known Drucker-Prager [1] and Mohr—Coulomb [2] models. These models can have either a static or a dynamic yield surface and an associative or non-associative flow rule. Multi-yield function models, as its name implies, consists of two or more yield surfaces, each of which can be static or dynamic. Modified Drucker-Prager/Cap model [28], having a static yield line and a dynamic cap, falls into this group. Initially used to model soil, this model is now widely used to capture hardening behavior of granular materials. There are two important innovations in this model compared to earlier models. First, the idea of a work-hardening cap is introduced as an alternative to the perfectly plastic yield surface such as the Drucker-Prager type or the Coulomb type of yield criteria. Second, the material relative density (or void ratio) is used as a state variable (or the strain-hardening parameter) to determine the successive loading of cap surface [29]. A number of modifications to extend this model have been attempted, but they all basically accomplish the same purposes as follows. 16 1) The dynamic cap bounds the yield surface in hydrostatic compression, providing the necessary inelastic hardening mechanism to represent the evolving plastic nature of the compact. In other words, volumetric plastic compaction when yielding on the cap causes hardening (under triaxial stress). 2) The static shear line controls the volume dilatation when the material yields in shear by providing a softening as a function of inelastic volume increase created as the material yields on the shear failure (and transition yield) surface. Unlike volumetric plastic compaction, volumetric plastic dilatation when yielding on shear failure surface causes shear-induced softening. This model captures some additional key characteristics of the powder, separating the influences of triaxial stresses and those of shear stresses during the compaction process. This model is adopted and modified in view of its simplicity, robustness and abilities to model ceramic powder essential material behavior. More details on the model and real material behavior will be presented in Chapter 2. In the review by Aydin et al. [9], other models of plasticity include Johnson et al.’s two-intersecting yield surfaces, Venneer’s double-hardening model, the well- known Cam Clay model and the semi-continuum semi-discrete model of Gurson l7 and Tvergaard. Some of these models are quite complex, while others lack the robustness of incorporating experimentally observed parameters. Many researchers have approached the problem of compaction using the plasticity model. Their efforts are described in details in Section 2.1. One particular research group, namely Mori et al. [30-33], is worth mentioning at this point. Mori et al. employ a rigid-plastic finite element method to simulate the compaction and sintering process of ceramic powder compact to predict the shape of the sintered product. The yield criterion for porous materials proposed by Shima and Oyane [17] is adopted to account for density changes. This yield criterion is originally based on the observation of metal powder and therefore may not be appropriate for ceramic powder. Mori et al. [30-33] also use a visco- plasticity constitutive model that relates stress g_ = CeéP with plastic strain rate as opposed to the conventional elastic-plastic constitutive g‘ = Ce 3" . However, their use of an iterative scheme to modify the geometry of the compaction die in the simulation for design purposes is novel. One of the main problems in the theory of plasticity is to determine the nature of the subsequent yield surfaces. This post-yielding response is described by the hardening rule, which specifies the rule of the evolution of the loading surfaces during the course of plastic deformation. Indeed, assumptions made concerning the hardening rule can introduce a major distinction among various plasticity models. 18 In determining the correct yield surface, Shima and Mimura [34] proposed two yield criteria based on their observation in their silicon nitride ceramic powder compacting experiments. Unlike the standard triaxial tests normally employed to determine soil mechanical properties, they have adopted a novel experimental procedures that enable three-dimensional independent loadings. Kim et al. [9] and Bortzmeyer [6] both conducted the standard triaxial tests in determining the mechanical properties of zirconia powder. Aydin et al. [7-8] combined data from a simple cold-die compaction with assumptions based on empirical approximations of powder behavior to deduce the characteristics constants of their alumina powder. 1.3.2 NUMERICAL METHODS Various numerical modeling methods such as discrete-element method (DEM), finite-difference method (FDM), boundary element method (BEM) and finite element method (FEM) may be applied for powder compaction modeling, as summarized by Aydin et al. [20]. The discrete element approach has the microscopic appeal of modeling individual particle interaction with its neighboring particles. This approach considers the powder as a discrete assemblage of particles with finite sizes and predicts the mechanical behavior of powder using governing law deduced from basic physical laws. However, DEM is extremely 19 computational costly and may not be feasible to simulate the behavior of more than a few thousand particles, considering that a powder compact of current concern consists of millions of particles. Though DEM provides insight into the behavior of the particle assembly, additional understanding of the particle-particle interaction during powder compaction at microscopic and macroscopic levels is required to fully establish the model. These numerical challenges posed by DEM are considerably less significant in all the other continuum approaches, such as BEM, FDM and FEM. In these approaches, the powder system is reduced to manageable sizes and treated as a continuum. FDM loses its appeal when it comes to problems with material nonlinearity, time-dependent material properties and complex arbitrary shapes. Similarly, BEM is considerable problematic in solving nonlinear analyses with a large surface-to—volume ratio, where the displacement and stresses at many interior points are desired. On the contrary, FEM is cost-effective in investigating details of powder processes that have direct impact on the quality of the compact such as density distribution, dimensional variations and mechanical properties. Furthermore, the obtained data can be extended for subsequent sintering analyses. 20 1.4 LITERATURE REVIEW OF SINTERING MODELING 1.4.1 SINTERING PHYSICS To model the process of sintering, it is necessary to first understand its kinetics. According to Barsown [11], the macroscopic driving force in sintering is the reduction of excess energy associated with surfaces by: a) the increase in the average size of particles (coarsening) and; b) the elimination of solid/vapor interfaces and the creation of grain boundary areas, followed by grain growth (densification). Coarsening and densifications are usually in competition. To obtain near- theoretical densities, coarsening has to be suppressed until most of the shrinkage has occurred. This is very difficult and requires great care because the driving force for sintering is quite small, usually on the order of a few Joules per mole, compared to a few kilojoules per mole in the case of chemical reactions. When coarsening happens, the free energy for sintering is expended in forming large grains, which are kinetically very difficult to remove due to their thermodynamic stability. Sintering happens because of five basic atomic mass transport mechanisms: a) Evaporation-condensation (coarsening) b) Surface diffusion (coarsening) 0) Volume diffusion (densification dependent on location of the source) 21 d) Grain boundary diffusion (densification) e) Viscous or creep flow (densification) To sum it up, for densification to occur, the source of material has to be the grain boundary or region between powder particles, and the sink has to be the neck or the pore region, such as in mechanism c), cl), and e) above: 1.4.2 SINTERING KINETICS The topic of sintering is a broad one that a thorough investigation will require more work than the author has set out to accomplish. Despite the benefits of researching the sintering physics, our concern in the current modeling process will be solely on how a compact sinters under different thermal or compressive loads. Thus, it is assumed that all the necessary conditions for densification by sintering are met, and will therefore not be concerned with issues such as the causes for coarsening, or pore-grain boundary interactions, or even the effects of curvature on partial pressures and vacancy concentration, except acknowledging the fact that their differentials produce the driving force needed for sintering. However, the current modeling efforts will entail basic understanding of sintering kinetics that is dependent on many variables, listed and elaborated as follows [11]. 22 a) m 0) Temperature: Increasing temperature enhances sintering kinetics, when more diffusion mechanism are activated. Green density. The higher the green density, the less pore volume that has to be eliminated, decreasing the sintering (or densification) rate. Uniformity of green microstructure: Non-uniformity caused by agglomerate can lead to abnormal grain growth. Fortunately, this is not our concern because the problem of agglomeration can be avoided by choosing the right powder to start with. Non-uniformity after cold die compaction that will lead to shape distortion, however, is within the scope of the current investigation. Atmosphere: In some cases, the atmosphere can enhance or reduce diffusivity. Most of our sintering efforts involve conventional furnace, and thus the effect of atmosphere is not a concern at this point. Impurities: Sintering aid that forms a liquid phase or low-temperature eutectics can enhance sintering kinetics. They suppress coarsening by reducing evaporation rate and lowering surface diffusion, and even suppress grain growth and lower grain boundary mobility. Some will go into solution and create vacancies on the sub-lattice that should in principal enhance sintering kinetics. This is not a major concern since high purity powder is being used in this work. However, the use of impure powder may be incorporated into our model by adjusting material parameters in future efforts. 23 f) Size distribution: Narrow-size distribution will reduce or eliminate abnormal grain growth. However, broad size distribution enables denser initial particle packing [35] and thus reducing sintering rate. The effect of size distribution is also investigated by Ting et al. [36, 37] and Darcovich et al. [38, 39]. 9) Particle size: The finer the particles, the more surface area and thus in principal favors sintering. In practice, however, very fine particles may lead to agglomeration. The agglomerates tend to sinter together into larger particles, which not only dissipates the driving force for densification but also creates large pores between the partially sintered agglomerates. Only the effects of particle size on the sintering rate and not its adverse effects of introducing agglomerates will be considered since abnormal grain growth is not present in our sintered compacts. Barsown [11] reported the work of Coble, who proposed 3 stages in sintering, during which different sintering kinetics take place: a) At the initial stage, the inter-particle contact area increases by neck growth, with relative density increases from 60 to 65%. b) At the intermediate stage, continuous pore channels are seen at the triple points (three-grain junctions) and relative density increases from 65 to 90%. c) Final stage: Pore phase is eventually pinched off. Continuous pore channels disappear. Individual pores are either Ienticular on grain 24 boundaries, or rounded within a grain. Pore and grain boundary mobility increases dramatically at this stage, and has to be controlled to achieve theoretical density. Many other researchers have modeled the process of sintering in like stages in their works thereon. 1.4.3 SINTERING MODELS 1 .4.3.1 BACKGROUND Preliminary qualitative models [11] proposed for the initial stage include lattice diffusion model and grain boundary diffusion model. These models relate the ratio of neck-to-sphere radius to the atomic volume, temperature, diffusivity, and sintering time. These models capture the experimental observation that grain boundary diffusion is preferred over lattice diffusion for smaller particles while lattice diffusion is favored at long sintering time, high sintering temperature and larger particles. Most densification, however, happens at the intermediate stage. This stage is strongly dependent on the details of packing, and therefore most difficult to model. Barsown [11] referred to Coble’s work, who by making simplifying 25 assumptions on packing geometry and densification kinetics, proposed relationship between fractional porosity to sintering variables for both volume diffusion and grain boundary diffusion models. These flux equations can be integrated to yield the rate of geometry change. Caution has to be taken when applying these models due to its many broad assumptions, such as the ever existing driving force to shrink the pores, or the proportionality of the average diameter of the sintering particles with a constant length in an assumed packing geometry. At the final stage, densification slows down significantly, and grain growth becomes significant. Because of the difference in partial pressure, grains that are smaller than the average will evaporate away, while those that are larger will grow with time. Grain growth models will be discussed in Section 1.4.4. Swinkels and Ashby [40] questioned the validity of the Coble’s immediate stage approach in the case of sintering of spherical particles. In their work, they omit the intermediate stage because pores in an aggregate of sintering spheres are not cylindrical. In the transition from initial stage to final stage, the pores form connected channels, but the channels vary in section and surface curvature, making the cylindrical approximation by Coble a very poor one. They link the initial stage to the final stage with a transition region by interpolating between them to make them continuous. 26 1 .4.3.2 POROUS VISCOSITY MODEL Based on a rigid-plastic framework, Mori et al. [30-33] models the sintering rate as a sum of normal shrinkage that is a function of relative density, and the plastic deformation rate induced by the difference in shrinkage strain caused by non- uniform density distribution in the compact. They also try to predict the sintering fracture initiated by tensile hydrostatic stress. Lippmann and Iankov [41] perform a numerical investigation of powder metallurgy forming processes based on a rigid-plastic finite element model. They also propose a coupled sinter-compaction model that includes the effects of local capillary forces. Unlike Mori et al. [30-33] or Lippmann and Iankov [41], most other researchers [42-50] have attempted the problem of sintering using the porous viscosity approach. In their two-part paper, Bordia and Scherer [42-43] justify the use of pure viscosity theory in developing the constitutive equations for porous ceramic compact. They compare various models [44-48] for the shear viscosity and bulk viscosity that are available in open literature in the framework of pure viscosity. All these models described Viscosities in term of porosity (or relative density) and sometimes grain size, and differ only in the nature of the dependencies. Bordia and Scherer conclude that the viscous Poisson’s ratio of the correct model has to be positive, which conforms to experimental observation. They show how a negative Poisson’s ratio will result in erroneous prediction of sinter shrinkage and fictitious high stresses. 27 Porous viscosity approach has gained popularity among researchers because of its simplicity and capability. Based on this approach, Riedel and Sun [49] develop a simple model in which a two-dimensional hexagonal grain structure is considered with pores at the triple points. The model assumes different pore sizes at the triple points by the introduction of a distribution function, but does not consider grain growth. The only transport mechanism considered is grain boundary diffusion. This model is later adopted by Zipse [10]. Du and Cocks [50, 51] also develop their own constitutive model by including grain growth rate explicitly in the porous viscosity law. On the other hand, Shinagawa [52], Kraft et al. [53, 54] and Kwon et al. [55, 56] extend the viscosity approach to model sintering by grain boundary diffusion while taking account of grain growth, density and temperature. Kim et al. [58-60] then include the power law creep mechanism in addition to grain boundary diffusion in their sintering model. Cocks [61], on the other hand, examines the theoretical aspect of the constitutive laws for the sintering in which power law creep and grain-boundary diffusion are the dominant mechanisms of deformation and densification. These models are useful for the understanding of the sintering physics, but may not be necessary for the effective control of the densification process. Kim et al. [62] thus propose an alternative approach that consists of only analytical expressions curve-fitted directly from the results of sintering experiments, and thus a constitutive equation describing precisely the real behavior of the material. 28 They [63] later implement the constitutive model in a finite element simulation to help design a near net shape processing of a relatively complex shaped alumina powder component. Thus far, the constitutive models for viscous sintering described above are mostly porosity dependent. A few researchers have adopted more novel approaches. Jagota et al. [64], for example, propose an isotropic model that depends on contact area. The constitutive relationship is mostly based on the behavior of individual inter-particle contacts, represented by contact Viscosities and the statistics of the packing. Unlike the other models, this model allows the description of low effective viscosity in the early stage of sintering and the difference between how coarsening and densification sintering mechanisms change the effective viscosity of the packing. However, the model lacks capability when the sintering body is almost fully densified, as several assumptions adopted in the model are no longer valid during that stage. On the other hand, Aydin et al. [65] develop a finite-element procedure to predict shape changes of the green powder compact based solely on the principle of mass conservation, disregarding any constitutive rules. Although the quantitative results may not meet strict industrial tolerance, the approach has demonstrated its effectiveness and qualitative predictive capability. 1 .4.3.3 SINTERING POTENTIAL 29 Sintering is basically a reduction of excess energy associated with surfaces by the elimination of solid/vapor interfaces and the creation of grain boundary area, followed by grain growth. Sintering potential is, therefore, the free driving force for sintering due to the minimization of this interfacial energy of pores and grain boundanes. As noted in Cocks’ work [61], Ashby comes up with expressions to compute sintering potential or sintering stress as . He divides sintering into two stages, with the first stage similar to Coble’s (as reported in Barsown [11]) initial and intermediate stages, and the second stage equivalent to Coble’s final stage. Assuming grain size is two times the particle radius on fine powder, G=2R, as =673 02 M =-3£D2 M at stage1 (D<0.9) G l-Do R l—Do (1.1) 0‘ mfg at stage 2 (D>0.9) (1.2) as =3LL where r=|: r where G is the mean grain size, D is the relative density, ys is surface free energy per unit area, r is the pore radius. Sintering potential in (1 .1) and (1.2) are apparently functions of D and G. The two different functional relationships are to account for the different porosity structures at the two assumed stages: 30 1) An array of spherical particles of uniform size at stage 1. 2) Spherical pores become isolated at 4-grain junctions at stage 2. Kwon et al. [55] observed that if R < 1 micron, (1.1) and (1.2) predict too large a value. With rs =0.9 J/m2 for alumina, as =1 MPa is valid only when R is about 5 micron, if equations (1 .1) and (1.2) are used. According to Cocks [61], Ashby had not considered the contribution of grain boundary energy to the sintering potential. Ashby assumed that the sum of the grain-boundary and surface areas around a given grain remains constant. But as material densities, the total increase in the grain boundary area is half the decrease in surface area. Modifying (1 .1) and (1.2), 03:33—(2—K-‘fll02 ZD-DO atstage1 (1.3) G 1-D0 l a =— +8 6 atstae2 1.4 . G [60%) cos() 9 ( > n, is the grain boundary energy; G is the mean particle size per unit area; 6 is the dihedral angle and 17,,(6) is given by Raj and Ashby [66] for two-grain-, three- grain- and four-grain-junctions. 31 2-grain: Fv (6) = 2?” (2 -— 3cos 6 + cos3 6) (1.5a) 3-grain: l Fv (6) = 7r — 23in"l[lcsc 6) + -1--cos2 6(4 sin2 6 — 1F — cos"l coti9 cos 6(3 - cos2 0) 2 3 J3 (1 .5b) 4-grain: _ 1 _ — a -— 2 Fv(t9) = 8 £-cos‘l 2 €086. (3 A} 3 Amnfl l A2 + Acos (43in2 6i — A2}? -— - 4cosl9(3 —cos2 0)sin-l( . J J5 Zmna 2 l where A = 3[s/2(4 sin 2 I9 — l)2 - cos 0] (1.50) On the other hand, Shinagawa [49] proposed a simple expression 1 as = D” 27 D (I - 0°)3 for both stages (1.6) {R Do (I‘D) where the material constants are 4 =0.5 andN =5. Apparently, these material constants are assumed to account for the geometrical details in powder packing. Because of the arbitrariness of the available models, Du and Cocks [50-51], and Brook [as referred to in 50-51] have used 1 MPa for sintering potential regardless. The value of 1 MPa is consistent with the experimental findings of 32 Rahaman et al. [47] and Venkatachari and Raj [48]. It is worth pointing out that the alumina powder used by Du and Cocks [50-51] and Brook [referred to in 50- 51] had a mean initial grain size of Go = 0.3 pm, which is in the range of our experimental interest. We in the current work have assumed constant sintering potential. 1.4.4 GRAIN GROWTH MODELS Grain growth occurs when the average grain size increases with time as larger grains consume the smaller grains. These happen with the migration of grain boundaries, or atoms from a thermodynamically unstable region to a more stable region. Grains with more than six sides will tend to grow, while those with less than six sides will tend to shrink. Straight grain boundaries would be stable and stay straight, whereas curved b0undaries would be unstable and move [11]. A thorough investigation of grain growth is beyond the scope of this study. For the purpose of this thesis, it is sufficient to know the following [11]: . The presence of second phase (such as pores, inclusions, or solutes) in grain boundary can have dramatic effects on grain boundary mobility, and may be rate-limiting. Pores cannot enhance boundary mobility, and they may or may not reduce it. . As the pores shrink, the mobility of the boundaries will increase. 0 Pores do not always shrink. They can coarsen. 33 o The mobility of the grain boundaries may become so large that the pores can no longer keep up with them. The boundaries simply move too fast for the pores to follow and consequently unpin themselves. Reducing grain boundary mobility and/or increasing pore mobility can prevent this pore breakaway. This study will ignore the coarsening of pores or pore breakaway, since our ceramic specimens do not suffer from these effects. Ostwald ripening [11], a classical grain growth model that assumes only the intrinsic grain growth kinetics (and thus excluding the possible presence of second phase in grain boundaries) features a parabolic dependence of grain size on time. Various researchers have investigated grain growth and come to slightly different conclusions. Du and Cocks [50] assume grain size to be a function of relative density. Rahaman et al. [47], however, find that grain growth is independent of the relative density, but dependent on the sintering time, sintering temperature and the grain size. Furthermore, Venkatachari & Raj [48] find that grain growth depends not only on time but also strain, with higher strain rates leading to a higher rate of grain growth. They conclude that grain growth rate can be enhanced by deformation. Based on those studies, Kim et al. [58-60] propose that grain growth during high temperature deformation can be summarized asG =G's +Gd, where G, is the static grain growth under stress-free condition such as pressureless sintering and 34 Cd is the dynamic grain growth due to deformation during various high temperature processing. Since this thesis is concerned only with pressureless sintering in which no external loadings are applied, the effect of Gd may be ignored. This model captures the effects of temperature, current grain size, and the grain boundary diffusion mechanism in material constants, as explained in more details in Section 4.1.4. Though this approach is more mechanistic than physical, it is simple and can be curve-fitted to different experimental results. 1 .5 THESIS OVERVIEW In the following chapters, the compaction and sintering models adopted in this thesis are presented. The numerical implementations of the compaction model is detailed in Chapter 2, followed by finite element results and discussion of the compaction simulations in Chapter 3. In Chapter 4, the numerical implementations of the sintering model are revealed. The obtained finite element results are examined in Chapter 5. Chapter 6 concludes the findings and provides a glimpse into the future works. 35 Chapter 2 COMPUTATIONAL MODEL AND NUMERICAL SCHEME FOR COMPACTION SIMULATIONS 2.1 CONSTITUTIVE MODEL An elastic-plastic mathematical model is used to model the monotonic cold compaction of ceramic powder. The behavior of ceramic powder loose particles is different from that of the ceramic compact. Ceramic powder particles are inherently brittle and do not deform plastically. On the other hand, ceramic compact can be assumed to be a continuum that behaves elastically until it reaches its yield strength, upon which plastic deformation occurs. In the current development of constitutive equations, Cauchy stress and logarithmic strain measures are assumed. By strain rate decomposition d§=d§e+d§p, (2.1) dg is a differential change in the total strain, dge is a differential change in the elastic strain and dg" is a differential change in the plastic strain. In equation (2.1) as well as the rest of the text, the underscore will be used to denote a tensor quantity. Assuming small elastic strains, BW (2.2) 35 2:: 36 g_ is the stress tensor and W is the elastic strain energy potential. For the case of linear elasticity, g=£e£e=Ce(§_EP) (2.3) Q" is the forth order elasticity tensor and g is the total strain and g” the plastic strain. This is the basis of our constitutive model. 2.1 .1 MOTIVATION It is worth pointing out some of the key elements, or rather the differences, in the current study and in those of other researchers. These differences constitute the main motivation for the study. Of particular importance are the availability of material data, and the consistencies of these data in the open literature. Bouvard et al. [67], for example, did not model the compaction behavior of alumina powder from a relative density value low enough to be within the range of our experimental interest. Their experimental relative density data are in the range of 0.62 to 0.7, which is not directly applicable to our cold die pressing analysis whose relative density value lies between 0.3 and 0.5. The same can be said about the experimental alumina data found in Zeuch et al. [68]. In addition, there were many differences in the assumptions over the compact material behavior among various researchers. For instance, Kim et al. [9] account for the influence of the relative density on the powder material 37 mechanical properties, whereas Zipse [10], Westerheid et al. [69] and Aydin et al. [7-8] do not. Poisson’s ratio is assumed constant by [6, 8, 10], and a function of relative density by [9, 58, 60, 67]. Aydin et al. [7] applied many simplifying assumptions when deriving their material parameters from the uniaxial die pressing experiment. Mori et al. [30-33] adopted a yield criterion developed for metal powder possibly not suitable for ceramic powder. Kim et al. [9] choose to model the deformation behavior of the ceramic powder to resemble the critical state model, which is non-dilatational on the critical shear yield line, while many researchers [8, 10, 53, 67] have favored the modified Drucker-Prager/cap model that allows dilation on the shear yield line. One interesting finding by Bortzmeyer [6] is that under hydrostatic loadings, the deformation is not purely isostatic. In fact, the axial deformation is found to be always smaller than the radial deformation due to the anisotropy of the powder at the onset of compaction. A triaxial test without initial pre-compacting the powder by Shima and Mimura [34] reveals that under certain strain-controlled conditions, the measured stress can be hydrostatic while the deformation is rather deviatory. Both authors conclude that the observed behaviors are typical of highly compressible rigid granular materials. Barreling on an initially cylindrical sample upon purely hydrostatic compaction, as observed by Kim et al. [9], is a consequence of this behavior. These observations also suggest more investigation into the plastic potentials associated with shear yielding. 38 Inconsistencies in the findings among researchers also stem from the differences in experimental observations. For example, Bortzmeyer [6] conducted similar experiments to those by Kim et al. [9] but they reached different conclusions. Kim et al. [9] and many researchers have found or traditionally assumed the relative density dependency of ceramic powder compact elastic constants. However, according to Bortzmeyer [6], the Young’s modulus of his zirconia powder was found to be relatively independent of the relative density but dependent on an 1 equivalent stress defined as a = (p2 + q2)2. The reason for this discrepancy in observation is still unknown. Finally, since the measured density is uncorrected for the elastic deformation, Bortzmeyer [6] uses a modified yield criteria to detour the experimental difficulty in determining elastic parameters. He justifies this approach by pointing out a literature reference showing that a variation in the mean particle size causes only a shift in the hardening law of the caps (iso-density curves) but does not change their shape. Most other researchers [8-10, 58, 60, 67] do not mention the effort of correcting for the effect of elastic recovery on the measured relative density. In addition, Bortzmeyer's [6] experiment yields a straight cap for the material hardening yield instead of a ellipse or hyperbole as observed by many [9, 67-68]. This thesis addresses most of these issues by conducting our own modeling efforts. Much of the present work relies on material data in the open literature. 39 Simulation results from different modeling approaches have been our main guidance in choosing our modeling approaches. Some earlier efforts such as the first order semi-empirical differential slice model by Briscoe and Rough [70] predict experimental relative density distribution qualitatively but lack accuracy quantitatively. Aydin et al.’s [25] hypo-elastic model fails to show the characteristic high and low-density regions observed in experiment. Our findings indicate that the modified Drucker-Prager/cap model and a novel hyperbolic cap model seem to predict the compact behavior accurately on most ceramic powder materials. Details about the material and geometrical models currently used are detailed in this chapter. Numerical results and their discussions are presented in Chapter 3. 2.1.2 HYPERBOLIC CAP MODEL It is intuitive to first discuss the behavior of the ceramic compact under the general triaxial loading and shear loading, particularly within the scope of failure criteria. For the ceramic compaction modeling, the hyperbolic cap model originally proposed by Kim et al. [9] has been adopted. The modified Drucker- Prager/cap model is also used for comparison and verification purposes. 40 Experiments by Shima and Mimura [34] show that plastic deformation of a ceramic compact depends not only on the second stress invariant (as in Von Mises material such as traditional solid metals), but also on the first stress invariant. However, dependency on the third stress invariant is negligible [25, 34]. Bortmeyer [6] conducted similar experiments and reached the same conclusion. In the light of this observation, the following yield criteria have been adopted for our hyperbolic cap model. Following classical theory of plasticity, the stress states within the envelope of the yield surfaces are assumed to be elastic. Details of each criterion are now discussed. 2.1.2.1 LINEAR DRUCKER-PRAGER SHEAR LINE q A fo'(p,q) linear Drucker—Prager yield line ”i ., Figure 2.1. Yield locus of the linear Drucker-Prager function The first criterion, namely the linear Drucker-Prager shear failure criterion depicted in Figure 2.1, is employed to account for the material behavior under shear loadings. When the material stress state goes beyond the linear Ducker- Prager yield line, the compact is assumed to deform continuously without 41 changing its volume (non-dilatational), just like the material behavior at critical state in the critical state model [3]. The ceramic compact is assumed to deform perfectly plastically when this happens. Because of the non-dilatational deformation of the compact along this line, the flow rule has to be non- associative. Adopting an associative model along the critical state will lead to inaccurate prediction of dilation [3]. The linear Drucker—Prager shear failure yield criterion reflects some of the important characteristics of ceramic behavior such as elastic response at lower loads, pressure dependent yield strength, plastic flow upon yielding, and elastic recovery during unloading. Though there are apparent shortcomings in this yield criterion, as thoroughly discussed by Chen [4], this criterion with the compliment of cap yield criterion is a logical approach in attempting the compaction problem. The linear Drucker-Prager yield criterion has the form [21] f0'= —a'p+r—k'=0 (2.4) where p = --:15-g:_l_ is the hydrostatic stress. (2.5) 1_ is the second order identical tensor and q is the equivalent stress: I q:(_3.£:£)2’ (2.6) 2 where §_ is the deviatoric stress, such that g = g. + pl (27) 42 a', k are material parameters. 1 is the shear stress: ,zfizmzi (2.8) J2 is the second invariant of deviatoric stress tensor: 12 =%[(0'x -0'y)2 +(ay —az)2 +(c7Z -ax)2]+r§y +732, +ng (2.9) ex, 0),, 02,230,,ryz , sz are Cartesran stress components at a pornt In the ceramic compact. If a' is zero, (2.4) reduces to the Von Mises yield criterion. Since 2' is proportional to q, the yield function can be easily plotted in the p-q space (also known as the Meridional space) as in Figure 2.1. fl is known as material friction angle. (1 is related to the material cohesion or material constant k. The yield criterion (2.4) can be rewritten as . q . =—a +—:k 2.10 In the p-q space, fo=q-(tanfl)p-d=0 (2°11) Comparing (2.10) and (2.11), it is apparent that d = «[515 (2.12) tan ,3 = .552 (2.13) 43 2.1.2.2 CAP SURFACE The second criteria, a novel hyperbolic yield locus proposed by Kim et al. [9] and depicted in Figure 2.2, has been used to model the hardening of material under general triaxial stress loadings. The relative density of the material is chosen to be the hardening parameter for the material. This is the ratio of the current density of the compact to its theoretical density. The magnitude of the relative density physically represents the amount of void in the ceramic compact. q A . f1(prqu) hyperbOIIC ..........yield line fo(p.q) linear d$ Drucker-Prager yield line 1’ p Figure 2.2. Hyperbolic-cap and linear Drucker-Prager yield functions This hardening parameter is generally agreed upon by Shima and Mimura [34] and Bortzmeyer [6] following their observation on ceramic compact behavior under compression. In fact, Shima and Mimura [34] has proposed a yield criterion with relative density being the hardening parameter. Other researches such as Mori et al. [30-33], Bouvard et al. [67] and Kim et al. [9] have adopted this proposal in simulating compaction of ceramic compact. On the other hand, Zipse [10], Westerheid et al. [69] and Aydin et al. [7-8] sidestep the direct use of hardening parameter by employing the modified Drucker-Prager/cap model that relates cap expansion to the total volumetric plastic strain. The volumetric plastic strain, however, can be related to the relative density of the compact. In employing the modified Drucker-Prager/cap model, the relationship between the hydrostatic pressure, p, and the total volumetric plastic strain is sufficient to represent the effect of hardening. Due to its shape in the stress-space, the hardening yield surface is generally known as the cap surface. The cap yield surface basically provides a measure to account for the hardening of ceramic compact under triaxial compaction. The cap surface moves out as the material yields under compression, thus increasing the yield strength of the material on subsequent loading. To capture the hardening of the powder compact under general triaxial compression, the hyperbolic-cap yield criterion of the form proposed by Kim et al. [9] is first adopted. fl (p,q,D) = 0 = i6 + cosh{A(D) . L6} — 3(1)), (2.14 a) 10 10 where A(D) and 8(0) are empirical values obtained from triaxial compression experiments by plotting the iso-relative-density curves at various relative densities in a graph of effective stress q versus hydrostatic press p. D in equation (2.14 a) is the relative density. To obtain a ceramic compact at its theoretical density under cold compaction is near impossible because ceramic powder Particle, being inherently brittle, doesn’t deform plastically and its shape is near 45 spherical. Previous experiments have shown that the relative density of a fully compacted ceramic preform in green state usually lies in the range 0.4 to 0.5 of theoretical density. The use of hyperbolic cosine function here provides a limit for D and is consistent with the physical reality. Care has to be taken with the use of hyperbolic cosine function. Since p and q in the current analysis are expressed in unit Pascal, the normalization factor 10° has been added to avoid exceeding the numerical limit of the function. Note that A(D) and B(D) is dependent on the choice of the unit for p and q, and therefore not dimensionless. It is obvious from (2.14a) that the hyperbolic-cap yield criterion does not have a unique solution at q = 0, or under purely hydrostatic loadings, since the differentiation of the yield function with respect to stress tensor is not smooth across the p-axis. This translates to plastic strains not having physically unique directions. However, this non-smoothness is what Shima and Mimura [34] observed in their experiment compacting the silicon nitride powder. Some researchers [9, 71] have attempted to resolve this shortcoming by introducing an artificial circular-arc yield criterion at an arbitrary low q value. This work will avoid this arbitrariness, which is perhaps unnecessary owing to the fact that the powder is under monotonic compressive loading in the concerned problem, which can be accounted in a numeric scheme. As such, the yield curve direction is unique as approached from the compressive region (or positive q-plane). 46 2.1.3 MODIFIED DRUCKER-PRAGER/CAP MODEL For verifications, the author’s have employed modified Drucker-Prager/cap model available in the internal library of ABAQUS 6.3 commercial codes to simulate ceramic powder compaction. Details of the model can be found in ABAQUS Theory Manual [28]. Figure 2.3. Schematic of flow for the modified Drucker-Prager/cap model in p-q space. The yield surface for the Drucker-Prager/cap model is represented in Figure 2.3 [28]. It consists of F, (the shear failure (yield) line), Fe (the cap yield curve) and I F, (a transitional yield curve for smoothening F3 between Fe). Similar to (2.11), the shear yield line is Fs=q-(tan,B)p—d=0 (2.141)) 47 where d is the material cohesion, and ,6 is the material friction angle. The plastic flow along the shear failure line is non-associative. Just as in the hyperbolic cap model, plastic hardening occurs when the powder compact yields on the cap. As the material hardens, the cap curve expands gradually from the dashed curves to the solid curve as illustrated with a dashed arrow line in Figure 2.3. The shape of the cap is controlled by a shape factor R (03, since the cohesion d is defined as d =[1-%tanfl)ac[28], where o". is the uniaxial compression yield stress. This is not seen as a limitation since it is unlikely this will be the case for real material. 55 When the material is non-dilatational, u/ = 0. Setting (11:19 will reduce (2.23) to the original model proposed by Drucker and Prager [1]. Following (2.22) and w = 0, the potential flow is assumed so that gP=EP§£=Epfl=Epg (2.24) 89; 8g where E” is the equivalent plastic strain rate defined as 01 7" = germ-P (2.25) i and ti: _ 2g g (2.26) 2.2.2.2 YIELDING ON MODIFIED HYPERBOLIC CAP As for the cap yield surface, it is assumed that the normality condition is fulfilled such that the flow is associative. In this case, yield criteria will the same as the flow potential. Furthermore, the flow potential is assumed to be such that it can be expressed in terms of 1) hardening parameter D and 2) either the stress tensor or scalar quantities p and q. The flow rule is then [72] 56 33 3f 18f Bf d P=d2—=dA—=d/t -——I — g 39. 32 I 3ap'+aqfl] (2.28) The plasticity model is completed by describing the evolution of parameter D in (4) with continuing plastic straining [72] d0 = h'(d§_p._0;.D) = h(d§_p.p.q. D) (2.29) For rate independent material, h must be homogenous of degree one in dg”. The plastic strain directions as determined from associative flow rule for the hyperbolic yield curve can be found in Appendix A. The following development has been adopted from the work of Aravas [72]. Subscript t indicates the beginning of an increment whereas subscript t+At indicates the end of increment. In an incremental approach, (2.28) can be rewritten as 1 Ag” =§A5pl+A5qflt+m where A5,, = —M[2f—J and Aeq = 154%) 3P t+At 34 ”A. Using equations (2.16), (2.19), (2.20), (2.30) and (2.31), ”HA: = 0e " KAgpl " Zmb'qnwm h _ 3 w are gtm _— £ =—£ A; zqe e zqt-l-At H 57 (2.30) (2.31 a, b) (2.32) (2.33) Equation (2.33) is a consequence of the flow rule being independent of the third stress invariant so that 524%, and g, +13: are coaxial, and linear isotropic elasticity assumption in (2.16) such that gem”, and g , +13, are coaxial. BUI (Tn-A: = —pt+At l +§qz+m Him: (234) 1 Comparing (2.32) and (2.34), Pr+At = pe + KAEP (2.35) and qM, =q, —3mgq (2.36) Finally, equation (2.29) in an incremental sense is AD = h. (A§p ’ 9'4“]; 9 Dt+At ) = h(A£p 1 [leg . Pr+At . 9r+At 1 Dt+Ar) (2-37) The problem of integrating the elastic-plastic equations then reduces to the solution of (2.35), (2.36), (2.37) and the following basic equations (2.38) and (2.39) _ 3f) [8f] _ Y — A8 — + A8 — — 0 (2.38) l {3‘1 t+At q 31’ r+At Y2 = f = f (17.20) = 0 (2-39) Equation (2.38) can be obtained by eliminating M in equation (2.31). The basic equations have to be fulfilled to correctly implement the constitutive model. Aravas solves equations (2.35)-(2.39) using Newton’s method. As such, the basic equations (2.38) and (2.39) for Newton’s iterative method become linearized equations [74] 58 dY, = —Y, (2.40) and de = -Y2 (2.41) Equations (2.35)-(2.37) and (2.40), (2.41) can then be grouped together to form A1 laAsp + 11123qu = b1 AZIBAEP + AzzaAeq =b2 (2.42 a, b) The coefficients A11,A12,A21,A22,b1, b2 are derived in Appendix B. Using this Newton’s iterative method, Asp and A5,, are set as zeroes at the beginning of the first iteration. At the end of each iteration, A5,, and A2}, are then updated as follows: Asp —> Asp +8A£p A84 —9 Aeq + 8A£q (2.43 a, b) The stresses and hardening parameter D (relative density in our case) are updated accordingly using equations (2.42)-(2.43). The iteration is performed until A5,, and A8,, both converge. Applying the principal of conservation of mass and assuming small strain analysis, the governing equation (2.37) can be reduced to (2.44) D = —Dé,, = h{£‘p,é‘q,p,q,D} (2.44 a, b) AD = —DAs,, = h{A£p ,Ae‘q , p,q, D} Thus 0”,, = D, + AD = D, — DMAek’; 59 D... = D' (2.45) 1+ A8; The implicit assumption in equation (2.44) is the relative small volumetric elastic strain such that it does not produce a significant change in the relative density of the compact. This is found to be the case in ceramic powder [7]. Irreversible plastic deformation, on the other hand, densities the compact under compression. This differs from the traditional plasticity assumption that volume change cannot be caused by plastic deformation. 2.2.3 LINEARIZATION MODULI In implicit code, finite element solver requires the specification of the linearization moduli at the end of the each increment for the computation of strain tensors tor the next increment, while fulfilling the equilibrium and compatibility equations. The linearization moduli is different from the material stiffness matrix in elastic- plastic model in that W 9". eé’ (2-46) II In 9'. = ”2' (2-47) IO 60 where the upper dot denotes rate quantities such that g, g" and g are stress rate, elastic strain rate and total strain rate. Q8 is the stiffness moduli as in (2.16) whereas gap is the linearization moduli. 2.2.3.1 LINEARIZATION MODULI FOR SHEAR YIELD SURFACE The linearization moduli for the non-associative linear Drucker-Prager yield surface assumes that of a traditional Von Mises yield criterion for g =qwhen IV = 0 is assumed. Thus [28]. 3.. A6")!- = Mam“, +2p—g—Aéij dull-ELM), (2.48) q. q. q q or more generally in terms of Lame constants, astound in Simo and Hughes [16] ® 932 = KL®£+2#(l_i)(l—l_l_®l)—2fl(l—i) 3_§e__2_§_e (2.49) 1 is the forth order identity tensor. The definitions of all the symbols used and the derivation of the linearization moduli itself are consistent with the current text and can be found in ABAQUS Theory Manual [28]. The derivations upon modification for the current application are shown in Appendix C1. 61 2.2.3.2 LINEARIZATION MODULI FOR HYPERBOLIC CAP SURFACE Govindarajan [71] derived the linearization moduli for the pressure-sensitive cap yield function in a way different from that of Aravas [72]. His work is summarized as below. Rewriting (2.20) using (2.31), the elasticity equation becomes Ut+At =£e[§r+Ar *Qrp _%A£p£_A£qflt+At) (2.50) and so a — e l anti-At a OHM _g 35%, ~38A£p1—8A£qgt+m -— A8,, 85—— 5M, —t+At (2.51) or m=ce[1_lfifl.1_ifln N_A£ ant-FA!) (252) ._ _ _ —t+ o a§t+At 3 a§:+A: aé’nA: q a§t+At It is shown in Appendix CZ that aEt+At 2.“ 3 1 567—; = —q_ ‘2']. T 2'11 T flr+At flr+At (2'53) —t+ e Using equations (8.1), (B2) and (8.3) in Appendix B as the basic equations again, the following set of equations results in solving for 3A5, and aAsq: C1131“ p + €12?)qu =(Dul + DIZ"t+Ar)'a§t+At (2-54) CZIBAEP + C228A£q = (0211 + Dzznt+m)' 85%, (2.55) 62 That is, [3MP = l [ C22 *C12][011£°3£1+A: + Dlzflt+At 'a§t+At] BAEq C11 C12 -C21 C11 0211332}er +D2221+At 'a§t+At C21 C22 (2.56) or in short 346 p = mprl'3§:+At T mpn EMA: 'a§r+At (2 57) aAEq mql _I_ ' a§t+At + mqn Hr+At ' a§t+Ar Since QHA, and 8;, +13, are known, equation (2.57) can be easily solved by grouping the coefficients C11,C12,C21,C22, 011101210211D221 as in Appendix C2. Equation (2.57), written in differential form is 3A6 p = mp1; + mpn EMA: (2'58) §r+Ar 8A8 3;”: = mg, l + mqn EMA! (2.59) Substituting (2.59) into (2.52), the linearization moduli is then ea)”, c e l 1 aflnAr .3;- = C P = C [l_§mplu—§mpn!.fl-mqlflH-All—manH-AIQHAI —A£ m] (2.60) 63 2.3 IMPLEMENTATION In the current implementation, material data has been obtained from different literatures. It is the intent in this work to experimentally obtain material data of the powder in use in the later stages of this research. Analyses are done on zirconia powder and its properties provided in the following sections. This material data has been adopted in simulations using the authors’ own code for the hyperbolic cap model and commercial codes for the modified Drucker- Prager/cap model available in the internal library of ABAQUS 6.3. 2.3.1 ADOPTED ZIRCONIA DATA Kim et al.’s [9] partially stabilized 3 mol percent Y203 zirconia powder (HSY-3.0) was purchased from Daiichi-Kigenso Kagaku Kogyo Co. Ltd., Japan. The chemical composition of the powder as provided by the manufacturer can be found in [9]. The mechanical properties of bulk zirconia they obtained from Ceramic Sources by American Ceramic Society are also repeated in Table 2.1. Table 2.1 Mechanical properties of bulk zirconia. Theoretical density (g/cm") 6.08 Young’s Modulus (GPa) 206 Poisson’s ratio 0.31 Kim et al. [9] conducted a triaxial compression test on the zirconia powder pre- compact and presented the relationship of stresses in the p-q space of the powder compact with its relative density on iso-density curves. They experimentally measured the Young’s modulus of their powder compact during the triaxial compaction test by measuring the variation of the axial strain with the axial stress in infinitesimal unloading. They believed that the Young’s Modulus is a function of relative density and expressed in equation (2.61). E = E0 exp[—(b(l — D) + c(1- D)2 )1 (2.61) where E0 is the Young’s Modulus of the bulk material, D is again the relative density of the compact. b and c are curve-tit parameters, where b = 12.6, and c = -6.99. They also used a mathematical model to estimate the Poisson’s ratio of the powder compact, as seen in equation (2.62). (2.62) 4 l+2(l-D)-3vo(l—D) v _1[4v0 +3(1-D)—7v0(1—D)] where v0 = 0.31. It turned out that the Poisson’s ratio is almost constant for the range of relative density in which compaction occurs (0.3= 0.42. However, the measured average relative density of the initial compact, or the cold die compact perform, is 0.31. Therefore, the yield criterion for relative density from this initial value D0 = 0.31 to D = 0.42 has to be extrapolated. These data have been adopted in the simulating zirconia powder compaction. As for verification through the modified Drucker-Prager/cap model, the following values have been adopted, also from Kim et al.’s work [9]. a = 0.03, ,6 = 54.3°, d = 1.5 MPa and R = 0.835 2.3.2 FRICTION BETWEEN POWDER BED AND DIE WALL Bouvard et al. [67] performed a sliding disk test between alumina powder and cemented carbide tool. They found that the Coulombic friction coefficient is approximately a constant of 0.18, and that the effects of compact density, its displacement and displacement rate on the coefficient of friction are weak. On the other hand, Zipse [10] conducted similar test by pressing green compacts of agglomerated alumina powder with 10% volume of Zirconia powder against tool steel rod without lubricant and found the coefficient of friction to be 0.23. They 66 therefore assumed constant friction in their compaction simulation. Aydin et al. [78] adopted a friction coefficient of 0.2 in [7] and 0.3 in [8] based on others’ research findings. Kim et al. [9] first obtained the experimental relationship between compaction and ejection pressures and then tried to find a match in their simulation by varying the simulation coefficient of friction. They found the best matching simulation friction coefficient to be 0.2. Briscoe and Rough [71] investigated the effects of zinc stearate lubrication on the friction between agglomerated alumina powder compact and the stainless steel die wall using a semi-empirical approach based on Jansen-Walker model [22] and the relevant experimental data. They concluded that the friction coefficient as a hyperbolic function of applied stress it the die is unlubricated, but is approximately constant at a value of 0.394 it lubricated. Considering the works of Aydin et al. [7] ( fl =0.2), Kim et al. [9] (,u=0.2), Zipse [10] (,u=0.23), and Bouvard et al. [67] (,u=0.18), a value of 0.2 has been used in this work for verification and simulation purposes. Carbon powder has also been used on the steel die as lubricant in this work. 2.3.3 FINITE ELEMENT APPROACH ABAQUS Standard 6.3 is used to simulate compaction analysis. The die wall, top punches and bottom punches are modeled as analytical rigid bodies. The following analyses have been performed: 67 1. In the preliminary simulation done for verification, the powder bed is modeled by a set of four-node axisymmetric continuum elements (CAX4 by ABAQUS reference), with the element stress order in the following sequence: a, , oz , 0,9, 1,2 , where the stresses are referred to a cylindrical coordinate system. The degree of freedom for axisymmetric elements are radial-displacement and axial-displacement. Compacts with two different aspect ratios are simulated, as listed in Table 2.2. Table 2.2 Geometry of simulated compacts. Height (mm) Radius (mm) Aspect ratio (H/2R) 4.114* 11.125 0.185 15.74“ 6.55 1 .202 *Geometry of specimen used in our lab ** Geometry of specimen used by Kim et al. [9] Both compacts are subjected to axial stress of almost 23 MPa, the loading applied in manufacturing our specimens. To compare with Kim et al.’s [9] experimental results, a load of 100 MPa is applied to the compact with higher aspect ratio. A 3-D geometrical is also built to simulate powder compact with embedded graphite pieces. For a general three-dimensional continuum 68 model, the stress order in an ABAQUS C3D8 element is ax, 0y, az,z'xy,z'xz,z' yr 3. To investigate the effects of embedded graphite piece on its vicinity, compaction of powder bed with the presence of graphite pieces in different orientations and locations are simulated using plain strain elements (CPE4) with stress order 0,, 0),, oz , rxy. This allows more detailed information than can be obtained using a SD computational costly analyses. 4. To investigate the differences between the cap model and the hyperbolic models, most of the analyses are run on both models. All the compaction simulations are broken into three steps. In the first step, the top punch moves down to compact the powder bed. In the second step, the compact is allowed to relax in the vertical direction while still in the die by removing the top punch from the compact. In the final step, the die wall is removed to simulate powder compact ejection. The final results represent that of a relaxed compact. The results of the analyses are presented in the Chapter 3. 69 Chapter 3 COMPACTION SIMULATION RESULTS 3.1 VERIFICATION Verification of the model is very important particularly when the material model is heavily dependent on experimental data. For this purpose, we adopt a compact geometry by Kim et al. [9], and compare their experimental results and our simulation results. These comparisons give us the confidence of the reliability of the model before proceeding to simulate other geometries of concern. 1 3.1.1 RESULTS AND DISCUSSIONS Kim et al. [9] perform uniaxial pressing experiment on zirconia powder, and through hardness tests, find the experimental relative density distribution as shown in Figure 3.1 (A). We have used'their compact geometry and applied the same axial loading on our quasi-static compaction simulation, as in Figure 3.1 (B) and Figure 3.1 (C). Since an axisymmetric analysis is involved, only the right half of the cylindrical compact is modeled. The dash line represents the axis of symmetry. Relative density distribution is indicated by a legend displaying the smallest to the largest values, marked with corresponding scale. All the 70 intermediate values are of equal intervals between the maximum and the minimum values. 2 I. R )1 I; 3 l# 4 l I \‘v I: I If i : 1 (k ‘1: : 5r" lid Density value .‘1 t I 3 9 I |( 7 I4 1: 0.405 I . I7 ‘ I 6\ 10: 0.460 ' I l T I I I l \2 W I—_\ \‘TT‘RK I I I7 ‘ 4 I I I-r/‘T‘TT/ l I L ' I {9" I f l I I /’ I ’5’ I 4:? . 4. g I”; I ( 13 I ( 112 r | ' 3.1 (A) 3.1 (B) 3.1 (C) Figure 3.1. Relative density distribution of the ejected right-half axisymmetric cylinder of height = 15.746mm and radius = 6.55mm from (A) Experimental results by Kim et al. [9], (B): the hyperbolic model (C): the modified Drucker- Prager/cap model. The disc is subjected to 100 MPa axial loading. It is fair to conclude that the experimental and simulation results, especially those of the hyperbolic model, are very much in agreement. The highest density region is at the top right corner, whereas the lowest density region is at the bottom right comer. This trend is consistent with the experimental results obtained by Aydin et al. [8,25], Briscoe and Rough [70], and Kim et al. [9]. 71 The differences in the simulation results between the modified Drucker- Prager/cap model and the hyperbolic model are primarily due to the different unloading mechanisms in the two models. The modified Drucker-Prager/cap model allows softening to occur, whereas the hyperbolic model assumes no softening when yielding under shear. These differences seem to lessen if the compact is subjected to a lower axial loading, as will be seen in the next section. 3.2 COMPACTION WITHOUT GRAPHITE PIECE Before introducing fugitive phase Into the powder compact, we first investigate the effects of aspect ratio by employing compact of different geometries. 3.2.1 RESULTS AND DISCUSSIONS The simulation results for the zirconia powder compact of two different aspect ratios, namely UR = 2.40 and UR = 0.37 respectively, are presented in Figure 3.2 and Figure 3.3. The results clearly show the differences in the relative density distribution between the compact with high aspect ratio and that of low aspect ratio, with 72 more uniformity observed in the compact with low aspect ratio. This perhaps explains the success of machining low aspect ratio disc to produce complex surface channels [75]. It has been known that given the same machining and sintering conditions, a compact with high aspect ratio has a higher probability to suffer distortion or even cracks. Density distribution non-uniformity in the compact with a higher aspect ratio will eventually lead to serious differential shrinkage and may result in product defects. The tool geometry and the friction between the die wall and the powder compact are the major reasons for the non- uniforrnity of the relative density distribution. Thus, lowering the tool aspect ratio and reducing this friction will help reduce chances of defects in the final products. 11 1 ((110 \ ([1, Relative density legend 2‘ I \9 \ \“‘~10 \- v49 1:0.3545 7:0.3814 7\ 3 1* \_ \ X). 8 \‘\ \\ “~~~ \ 7 \~.,\ 2 : 0.3622 3 : 0.3853 6\.\ \\ I‘m \A 6 \\ 3:0.3661 9:0.3391 4 : 0.3699 10:0.3930 '5 _ r” 4 5:0.3738 11:0.3968 2 . (51 {51 6.0.3776 3.2(A) 3.2(B) Figure 3.2. Relative density distribution of the ejected right-half axisymmetric cylinder of height = 15.746 mm and radius = 6.55 mm from (A): the modified Drucker-Prager/cap model and (B): the hyperbolic model. The disc is subjected to 23 MPa axial loading. 73 1’: I 3.3(A) 4 3.3(B) Figure 3.3. Relative density distribution of the ejected right-half axisymmetric disc of height = 4.114 mm and radius = 11.125 mm from (A): the modified Drucker-Prager/cap model and (B): the hyperbolic model. The cylinder is subjected to 23 MPa axial loading. Relative density value markers are the same as those in Figure 3.2. Von Mises stress (MPa) legend 1:0.009 2:0.158 3:0.307 4:0.455 5:0.604 620.753 7:0.902 8:1.050 9:1.199 10:1.348 34:3 2 Mg: \\ \1Il9 : 3 (>\ : l 3 3 4 [fix (‘1 46. ] fire: 3.4(A) ' 3.4(B) Figure 3.4. Residual Mises stress distribution of the ejected right-halt axisymmetric disc in Figure 3.3 (A): the modified Drucker-Prager/cap model and (B): the hyperbolic model. 74 As seen in the earlier section the modified Drucker-Prager/cap model and the hyperbolic model do not produce exactly the same density distributions in both cases as in Figure 3.2 and Figure 3.3, although with lower axial loadings the results seem to match better. Again, the highest density region is at the top right comer, whereas the lowest density region is at the bottom right corner. Figure. 3.4 again shows how the residual stress distribution differs slightly between the two models. The residual stress will play a significant role in sintering, affecting the sintering potential as explored in Chapter 4 and Chapter 5. The average density of the relaxed compact in Figure 3.3 was calculated using equation (15) 2WD,- 2": where i denotes each element. The simulation yields an average relative density (3.1) Dave = of 0.382 (hyperbolic) and 0.375 (modified Drucker-Prager/cap). whereas our measurement of real specimen yields 0.39, with a percent difference of about 2% and 4% respectively. 3.3 COMPACTION WITH GRAPHITE PIECE 75 To support the experimental interests of fabricating micro-channels, compaction simulations are done with graphite piece embedded into the powder bed. In our previous experiments, we have used a variety of the fugitive phase ranging from polymeric materials to graphite [76]. Polymeric materials are easy to shape into a desirable geometry, for example, using rapid prototyping. And we have had success in making the features such as channels and manifolds that are exposed to the surface of the powder compact. However, when they are embedded inside the powder to make enclosed channels and chambers, the stored energy during compaction is released during ejection to crack the compact. The stored energy is proportional to the size of the fugitive phase. It may be possible to make a smaller feature in microns and submicrons. Shin et al. [76] found success in making internal cavities using graphite as fugitive phase. Due to the high stiffness of the graphite, the stored energy is not high enough to demolish the compact. Another important aspect of graphite is clean- bumout above 600°C. Graphite of circular cross-section and rectangular cross section are used as the fugitive phase in the current study. Graphite piece of circular cross-section has not caused cracks in the samples. On the other hand, the use of rectangular cross-section has demanded more consideration due to the cracks that it may induce at the sharp edges. 76 Of particular interest to this work is the potential crack sites, which may well be predicted by our compaction models. Various researchers [6, 9, 25, 30-34, 53, 67-69] acknowledge how non-uniform density distribution leads to shape distortion and even cracks during sintering. But only a few [9, 25, 63, 67-69] have explicitly stated their view with regard to potential crack sites on green bodies or the green bodies themselves after compaction and before sintering. While some [25, 63, 67, 69] provide only a general comment on the poor quality of green compact for local variation of relative density, others [68, 9] specifically say that the cracks can form after die compression. Zeuch [68] hold the view that those cracks form at locations with high relative density gradients, whereas Kim et al [9] believe that they form at sites with high residual stresses. Finally, only Mori et al [32-33] have performed a scientific investigation on fractures after sintering on ceramic powder, and they developed a fracture criterion based on differential in shrinkage rates. And the fractures they study are those observable by naked eyes, not micro-cracks. 3.3.1 RESULTS AND DISCUSSIONS 3.3.1.1 THREE DIMENSIONAL COMPACTION 77 (a) Modified Drucker-Prager/cap model with zoomed in section Legend for both (a) and (b) 1: 0.3381 411: 0369656 . 4: 3’ g" 1“\\ 1| ‘ . {fl .1" 1‘ Figure 3.5. Relative density distribution after ejection for a quarterly powder bed with embedded graphite piece subjected to 11.5 MPa. The empty spot is where the graphite resides. Powder bed initial thickness is 6.0 mm, and its radius 11.125 mm. The graphite piece has a cross section of 1 mm X 1 mm. 78 Figure 3.5 shows the relative density distribution of a quarter model with graphite piece embedded in the powder bed. A full fine mesh of the three-dimensional model is very computationally costly. As a result, refined mesh in the region is used close to the graphite piece and a coarse mesh elsewhere. The empty spot is where the graphite piece lies. It is obvious from the figure that the highest density gradient occurs in the vicinity of the graphite piece, where cracks and micro-cracks occur in physical specimens. As a result, the information in this region will be most important in the current study. To save computational resources, the St. Venants’ principle has been applied to the three-dimensional problem, which reduces it to a two-dimensional one using plain strain analyses. 3.3.1 .2 TWO DIMENSIONAL COMPACTION An investigation is conducted on the effects of the embedded graphite piece on its surrounding powder flow. Specifically, we study how the size of the graphite piece, its orientation in the powder bed, and the loading applied to the compact affects the relative density distribution and introduces micro-cracks to the final sintered product. Studies are performed using the hyperbolic cap model and also the modified Drucker-Prager/cap model. 3.3.1.2.1 SIZE AND ORIENTATION EFFECTS 79 Simulation results from the plain strain analyses (Figure 3.6 (8)) indicate similar relative density distribution to the one in Figure 3.5 (as highlighted by the dotted box), but with greater details. The results for a loading of 23 MPa with graphite pieces of different sizes are presented in Figure 3.6 and Figure 3.7. Modified Drucker-Prager/cap Model _ Dvalue 1: 0.350 4’ , ~ 1: 0.331 12: 0.403 . . 120.397 12 \ , . 11 12, . \ j $97 . I I /2 1 n093 10937 (A) (B) (C) Figure 3.6. Relative density distribution after ejection simulated using the modified Drucker-Prager/cap model (top row) and the hyperbolic model (bottom 80 row). A graphite piece is embedded in different orientations and locations as in (A), (B) and (C). The graphite piece has a cross-section of 2 mm X 2 mm. All compacts are subjected to 23 MPa. The preform has a thickness of 6 mm. Modified Drucker-Prager/cap Model D value D value D value 1: 0.3271 1: 0.3237 1: 0.3699 11: 0.3570 _ 11:, 0.3550 . 1: 0.3453 11 A}; ,_‘ ' ’72: \\\ 9 1+8 I I I) :\.\' \ \ 9 (A) (B) D value D value D value 1: 0.3565 1: 0.3578 1: 0.3582 11: 0.3931 11: 0.3943 11: 0.3936 (A) (B) (C) Figure 3.7. Relative density distribution after ejection simulated using the modified Drucker-Prager/cap model (top row) and the hyperbolic model (bottom row). A graphite piece is embedded in different orientations and locations as in 81 (A), (B) and (C). The graphite piece has a cross-section of 1 mm X 1 mm. All compacts are subjected to 23 MPa. The preform has a thickness of 6 mm. The relative density distribution predicted by the modified Drucker-Prager/cap model and the hyperbolic cap model agree well at least qualitatively in the cases where the graphite pieces lie flat in the powder bed. In all cases, modified Drucker-Prager/cap model predicts lower relative density than the hyperbolic cap model. This is expected as the modified Drucker-Prager/cap model allows for plastic dilatation during unloading, and hyperbolic cap model does not. It is also observed in both models that high gradient in relative density occurs at the comers of the graphite pieces. Major deviations between the two models are more apparent in the case where the graphite pieces lie diagonally in the powder bed. These deviations are due, again, to the effect of unloading. The relative density distributions for both models after compaction but before top punch removal are very similar. During the top punch removal and compact ejection steps, the modified Drucker-Prager/cap model softening behavior starts to dominate, leading to the deviations shown in Figure 3.6 and Figure 3.7. The deviations are larger for the case of diagonal orientation, and show increasing trend with increasing compaction load, as shown in the next section. 82 3.3.1.2.2 LOAD EFFECTS Modified Drucker-Prager/cap Model 6\ \\ ;:-‘ t I "‘1 . \\ «3 ‘~ _ 1—D5U‘fl I \ .__ 4/ 1"- — . “A _ 7 L__ I -2. _— 4 7 l , 5 : 3x1 I "I 3 I 31/1’42 '[ 2 1 41' , " '1 1 I I -” / 2 1 I r ‘ /1 . I 4~ 1 1 A L ' a 1 (A) (B) (C) EH11 ‘3‘10 too: 0 ./-"5 '\ 7 7‘ l 4 . '3 ' J \ . ‘ I i ‘ 8 _ I\ f' ‘. ( 1 4 \ I. 6 / I \\ \K’ME .-.. y“ ‘x ‘k ‘ l t '. 9 5. D Value D Value D Value 1: 0.3338 1: 0.3347 1: 0.3351 11: 0.3653 11: 0.3642 11: 0.3637 (A) (B) (C) Figure 3.8. Relative density distribution after ejection simulated using the modified Drucker-Prager/cap model (top row) and the hyperbolic model (bottom 83 row). A graphite piece is embedded in different orientations and locations as in (A), (B) and (C). The graphite piece has a cross-section of 1 mm X 1 mm. All compacts are subjected to 11.5 MPa. The preform has a thickness of 6 mm. Observation of experimental results has shown that by reducing the compaction load, micro-crack formation can be reduced, or even eliminated. Thus, understanding the effect of compaction loading on the powder flow around the graphite piece is very important. The results of the plain strain analyses for a loading of 11.5 MPa (half of the previous load) with graphite pieces of 1 mm X 1 mm cross-sections are presented in Figure 3.8. Comparisons between Figure 3.7 and Figure 3.8 show that for both material models, the relative density gradient does reduce accordingly with the reduction in compaction load. 3.4 EXPERIMENTAL OBSERVATIONS Researchers [6, 9, 25, 30-34, 53, 67-69] have postulated that crack sites occur at regions with high relative density gradient. The presented modeling efforts have thus far provided us a qualitative assessment as to where potential cracks will form. For the model to be useful, a quantitative analysis has to be conducted to determine a critical gradient value that initiates the micro-cracks. This can only be accomplished with the aid of experimental findings. Efforts were taken to produce these specimens and compare them with simulations results. 84 3.4.1 SIZE AND ORIENTATION EFFECTS (A) (3) Figure 3.9. Micro-cracks observed in both samples with graphite piece of cross- section 2 mm X 2 mm subjected to 23 MPa compaction load. Graphite piece is oriented flat in (A) and diagonally in (B). Experimental observations in Figure 3.9 and Figure 3.10 reveal that reducing the size of the embedded graphite reduces the presence of micro-cracks. This can be explained by the fact that the stored energy in a smaller graphite piece is smaller, and thus the elastic recovery of the smaller graphite piece is proportionally smaller. A large enough elastic recovery can cause yielding at vicinity of the channel left behind by the burned graphite, leading to micro-crack formation. 85 Pencil leads of varying sizes have been successfully used to fabricate the channels using the same technique [76]. However, to fabricate a complex network of channels needed in our micro heat exchanger, the use of fugitive phase with more complicated geometry is essential. More importantly, experimental findings indicate the need to search for other alternatives that may eliminate micro-crack formation entirely. (A) (3) Figure 3.10. Micro-cracks observed in both samples with graphite piece of cross- section 1 mm X 1 mm subjected to 23 MPa compaction load. Graphite piece is oriented flat in (A) and diagonally in (B). 3.4.2 LOAD EFFECTS Simulation results presented above imply that reducing the compaction load can minimize micro-crack formation. Figure 3.11 shows us that by reducing the compaction load, micro-cracks in the sample with embedded graphite piece lying 86 flat are eliminated. In the sample with graphite piece lying diagonally, they are greatly reduced. Figure 3.11. Lowering the compaction pressure helps to reduce or even to eliminate micro-cracks in samples. 3.5 MICRO-CRACK SITE PREDICTIONS The relative density distribution is capable of predicting the location of potential crack sites of the samples by pinpointing where the gradients occur. However, the local gradient values computed at the critical corners have not been consistent enough to establish a critical value that aligns with our experimental findings, based on which crack initiation can be estimated. This difficulty in obtaining such critical value may be due to the inability to control the sharpness of the comers. Comparisons between the simulation results and our experimental findings above reveal the following: 87 1) 2) In the case where the graphite piece lies diagonally, the modified Drucker- Prager/cap model predicts a high relative density gradient around the corner encircled in Figure 3.7 (A), but not at the top or bottom comers. The hyperbolic model, on the other hand, predicts much lower gradients at the same regions. Experimental findings with different compaction loads indicate that at a low load (<1 klb), micro-cracks form at the top and the bottom comers. At greater loads (1 to 2 klb), micro-cracks will also form at the encircled region of Figure 3.7(A). In other words, if the relative density gradient is used as the sole criterion for micro-cracking, then the modified Drucker-Prager/cap model does not show predictions in line with the experimental observations. In the case where the graphite piece lies flat, both the modified Drucker- Prager/cap model and the hyperbolic model predict high relative density gradients around the comers. However, higher computed gradients are observed for the case of graphite piece lying in the mid-plane, compared to those shifted from the mid-plane. This contradicts our experimental observation because samples with graphite piece shifted from the mid- plane have consistently shown more tendencies to crack formation. In addition, the gradients computed at the corners of channels with the flat orientation are consistently larger than those at the top and bottom tip of the channel oriented diagonally. These gradients are erroneously lower compared to those computed by the modified Drucker-Prager/cap model at the mid corner of the diagonally-oriented channel. 88 Since we do not experimentally measure the local relative density distribution at the critical corners, we can account for these discrepancies only to one or more of the following: 1) Both models fail to predict the local relative density gradients with accuracy, even though the predicted global relative density distribution does help in locating potential crack sites. This failure may be caused by inadequate material modeling during unloading, unrealistic ejection boundary conditions or inadequate interface modeling between the graphite piece and the powder compact. 2) Relative density gradient may not be the best criterion, or the only criterion to identify the origin of micro-cracking. 3) There are uncontrollable parameters in the experimentations; the most common being the difficulty in aligning the embedded graphite piece perfectly. Having acknowledged the possibilities presented above, another route has been pursued to identify the origin of micro-cracking. Through the use of a flag in the user subroutine by which the hyperbolic cap model is implemented, it is observed that the sites of micro-cracking is very much coincident to the sites where shear yielding is predicted to be severe, often during the top punch removal. Since the modified Drucker-Prager/cap model is an in-built model in ABAQUS, efforts in extracting the same information can be technically difficult or even inaccurate as 89 similar flagging feature is not built-in. The returned stress states in the inbuilt model have to always lie within the yield surface. Thus, some guesswork may be necessary to obtain similar information from the modified Drucker-Prager/cap model. To understand the role of shear yielding, it is imperative to first understand its relation to the unloading of a sample compacted by uniaxial pressing. 3.5.1 UNLOADING UNIAXIALLY-PRESSED SAMPLE 02‘ q A A 82 C p Figure 3.12. Loading and unloading paths the unaxially-pressed sample, 2 being the axial direction. Aydin et al’s work [7] reports on ceramic compact behavior upon uniaxial pressing and unloading. According to Aydin et al [7], the loading and unloading 90 paths follow A-8 and BE respectively, as illustrated in Figure 3.12. ‘These paths assume constant elastic mechanical properties of the sample, which are not the assumptions of the current work. Nevertheless, it is sufficient for our purpose of understanding the effect of unloading. During unloading, the material stress states go from point B to point D in a straight line having a constant slope of K +§y in the oz —3z space, and a corresponding path BC in the p-q space having a constant slope of BL. «EX (Alternatively, path CD shares the same slope but with opposite sign since it is a mirror image about p-axis if BC is to continue below q=0. q>=0 by its definition in equation (2.6)). This is due to the assumption that elastic unloading occurs within the yield envelope, and that the powder material is isotropic linear elastic in the elastic region. When the stress state hits point D, shear yielding begins. As the hydrostatic stress is reduced, the stress state follows the non-linear path DE. Upon E, there is still radial stress within the compact (since the compact is not yet ejected), but with the top punch already removed, 0,, = 0. 3.5.2 SHEAR YIELDING AND CRACK SITES It is found in the current study that the compact material does not necessarily yield by shear during the top punch removal step, while locally most compact 91 material points will yield by shear after the ejection step. This is crucial since crack formation is likely to be related to shear yielding during the top punch removal. It is observed that the sites that start to yield by shear during top punch removal in the simulations and continue yielding through most of the ejection step are coincident with the sites of micro-cracking in our samples. Figure 3.13. Simulation results identifying regions that have yielded by shear (lighter area) before the ejection step under 11.5 MPa pressing with embedded graphite piece of the two sizes mentioned above. Figure 3.14. Simulation results identifying regions that have yielded by shear (lighter area) before the ejection step under 23 MPa pressing with embedded graphite piece of the two sizes mentioned above. Figure 3.13 illustrates the elements that yield by shear during the top punch removal step for the case of 11.5 MPa die-press load with embedded graphite pieces of the two different sizes mentioned. The lighter regions are the regions already yielded in shear before the ejection step. Note that the only case that has not yielded by shear before the ejection step is the one with smaller graphite piece lying flat in the mid-plane. This is also the case where our sample has no micro-cracks. For both the cases in which the graphite piece lies diagonally or flat in the mid-plane, the regions of shear yield in the simulation results are consistent with the origins of micro-cracks. Lastly, for the case in which graphite piece lies flat and shifted from the mid-plane, several samples have shown major crack detectable by naked eye near the surface of the samples. And we find that the simulation predicts this orientation to have yielded most during the top punch removal step. For the sake of completeness, we also include the simulation results for samples under compaction load of 23 MPa in Figure 3.14. Again, the negative effect of increasing compaction load is shown in the figure. 3.6 CONCLUSIONS The compaction simulation provides a valuable tool to understand the powder behavior during compaction especially with fugitive phase. The finite element simulations in this study have led to the following insights: 1) Machining ceramic compacts with a low aspect ratio was successful due to the uniformity in its density distribution. 2) The geometry and positioning of the fugitive phase can be designed by observing the regions of yield before ejection and the relative density distribution of the compact after ejection. 94 3) 4) 5) 6) Using a lower compaction load helps eliminate micro-crack formation. The relative density distributions predicted from the modified Drucker- Prager/cap model and the hyperbolic model help to identify potential crack sites. It provides useful visual information of the global distribution, and local sites where higher gradient occurs. However, the local gradient may not be accurate enough to be used to identify crack or to establish a critical value. Using a flag in the user subroutine umat, the hyperbolic model is able to identify regions where shear yielding occurs. Current results show that the regions of prolonged shear yielding (that is, regions that start yielding by the completion of the top punch removal step) most likely induce micro- cracks. These sites are consistent with the crack sites in real samples. Further investigation is needed considering the differences in the unloading mechanisms of the models and their effects on the simulation results. In particular, information regarding the shear yield and its relation with micro-cracks are essential in producing a better compaction model. 95 Chapter 4 CONSTITUTIVE MODEL AND NUMERICAL SCHEME FOR SINTERING SIMULATIONS 4.1 CONSTITUTIVE MODEL Due to the non-uniformity in the compact after unaxial pressing, different regions of the powder compact experience different sintering rates. This gives rise to strain gradients that include the volumetric and shear components. Densification happens due to the volumetric strain, while shear deformation relaxes the shear stresses in a compact. Thus the sintering behavior depends on the relative rates of shear and densification. A constitutive model is needed to account for the sintering behavior. Different sintering models are briefly described in Chapter 1. In this Chapter, the sintering model adopted and the reasons behind its selection are presented. 4.1.1 CONSTITUTIVE MODEL BACKGROUND Several mechanisms, such as diffusion, plastic flow, and power law creep contribute to densification. Since particle size plays a significant role in the 96 primary sintering mechanism in operation, it is important to first note the small particle size of zirconia powder used in this work as listed in Table 4.1 Table 4.1 Particle sizes of 3 mol percent yittria-stabilized zirconia Powder Brand Particle Sizes reported (micron) Cerac Z-2003 100 % < 3, 90% < 1.79, 50% < 1.23, 10% < 0.77 Tosoh TZ3YS 0.59 Nanbo SG3-TZP 0.2~0.4 Nanbo SG3S-TZP 0.2~0.4 According to Cocks [61], diffusion tends to be the densification mechanism in operative at low pressure and in fine grain ceramic materials. The primary path for diffusion can either be through the grain or along the grain boundary, with the overall rate of the process either controlled by diffusional rate or the rate of chemical reactions when material is added or removed from a grain boundary. Researchers such as Zipse [10], Kwon et al. [55-57], Kim et al. [58-60] have assumed that diffusion is the dominant material transport mechanism, which is true for a wide range of ceramic materials, particularly for oxides. In fact, in Kwon et al. [57] has shown that the effect of a non-diffusional mechanism, such as the power-law creep, is negligible when sintering happens at a low stress state as in pressureless sintering. Different densification maps provide further evidence on the dominance of diffusion [15, 57, 77]. Based on their findings, the powder listed in Table 4.1 can therefore be reasonably assumed to sinter under the diffusion in 97 our current investigation on pressureless sintering. Subsequently, our search for a constitutive model can be narrowed down to the one based on diffusion. It is important to note that the lattice diffusional flow known as Naborro-Herring creep and the grain boundary diffusional flow known as Coble creep are not the concern of this thesis. According to Helle et al. [77], these deformation mechanisms occur only when the grain size of particles is significantly smaller than the particle (i.e. smaller than the neck size). The powder used in this work is so fine that each particle is mono-grain. Shinagawa [52] summarizes with great precision the groundwork leading to the constitutive model employed here. He reports on Ashby et al.’s primitive two- sphere model in modeling the densification mechanisms of powder particles in HlPing (hot isostatic pressing), which leads to the development of densification rate equation. Shinagawa also follows up on McMeeking et al.’ [78] efforts on developing the general constitutive equations for the intermediate stage of sintering, taking into consideration the shear deformation in the two-sphere model. Artz [12] and Fischmeister et al. [13, 14] develop mathematically complex equations to describe monosized spherical powder packing, which are later simplified by Helle et al. [77] in correlating the mathematical model to their HlPing experiments. Du and Cocks [50-51] and Kwon et al. [55-57] have proposed constitutive equations for the densification and grain growth under diffusion. With 98 regard to the macroscopic behavior, Shinagawa [52] points out that metals and ceramics show linear viscous flow in pressureless sintering. 4.1.2 NEWTONIAN FLUID MODEL BASICS The Newtonian constitutive model has provided the framework to model sintering materials, including sintering glasses and some polycrystalline materials [42-43]. The model can be summarized below. Sij =2flpélj (4.1) tr(g) = 3 K p(tr(§') —36" f) (4.2) pp and K pare the shear and the bulk viscosity. Si} and éij are respectively the deviatoric stress and strain rate. The free strain rate 3, is the linear contraction rate of a freely sintering body (i.e. one that is not subjected to external load or other constraints). Equation (4.2) indicates that the volumetric strain rate is the 3K linear superposition of the free strain rate 33, and the stress induced rate I (f), r a tr(2) such that rr(§')=3éf +—3f_' This linear additivity implies the inherent p assumption in this model, namely the sintering mechanism is not affected by imposed stresses, even though the amount of strain rate is. 99 Rewriting (4.1) and (4.2) in terms of principal stresses and strains, where subscripts 1, 2, 3 denote the principal directions, results in . . l €1=€f+-E-_—[O'l—Vp(0'2+0'3)] (4.38) P . . 1 £2 =£f+E—[a'2—vp(al+a3)] (4.3b) P . . l 83 =€f +E—[0'3—Vp(0'1+0'2)] (4.3C) P E p and vp are the uniaxial viscosity and the Poisson’s ratio of the porous body. By analogy to Hooke’s Law, ,u = if— (4.4a) P 2(1+ VP) E K p = —£— (4.40) 3(1— 211p) 9K 1;” = _P_”P_ (4.4c) SKI, + up 3K — 2 — ( P ”P ) (4.40) vp - 2(3Kp +,up) The differences between the various constitutive models that have been proposed lie in the assumptions made about the dependence of the viscosity on porosity and grain size [42-43]. Theoretically, when porosity-e 0, the constitutive laws for incompressible fluid (included in Appendix D) should be recovered. 100 As porosity—> 0: 1) Ep -> 3,11,, and vp —> %, noted by comparing (D.3a-c) and (4.3a)(4.5a) 2) Kp—aooasv p —> g, as observed from (4.4b) (4.5b) E 3) pp —4 —3”—, noted by comparing (0.1) to (4.1) and (4.4a)) (4.5c) In principle, ,up and K p can have values ranging from 0 to co. However, sinter- forging experiments have always indicated vp > 0 [42-43]. Therefore, to conform with experimental observations 2;! K p > —3£ (4.6) A violation of (4.6) would result in an inaccurate or even erroneous modeling, as has been discussed in details in Bordia and Scherer [43]. One should note the response of model when the applied hydrostatic tensile stress is equal to sintering stress, such that 01 = 02 = a3 = as. In this case the sintering contraction stops as 5'1 = 32 = 33 =0. Consequently, averaging (4.3a) through (4.3c) yields . . . . as 5 =3 =6 =8 +—=0 4-7 1 2 3 f 3Kp ( ) 0 Therefore 3 = - S 4.8 f 3Kp ( ) 101 Equation (4.8) relates the sintering stress to the free sintering strain in a sintering material. 4.1.3 NEWTONIAN FLUID MODEL TO SINTERING MODEL In general, the problem of sintering can be summarized as follows. The total strain rate 8,-1- is decomposed into the elastic strain rate 35’ and the plastic strain . pl rate 31.]. . . l ., +5? (4.9) 11' The plastic strain rate 31.5.” is assumed to obey associated flow rule such that 3’” =éC'Q-l- 8's”! (4.10) (plu— g and _I_ are the second order tensors as defined in equations (2.26) and (2.5) respectively in Chapter 2. 8'" and é” are respectively the creep rate and the densification rate. Numerically, the stress strain relationship is 0t+At =Qd(€t+m "Epl —A£pl) (4.11) Equation (4.11) is similar in concept to the one seen in equation (2.20) of Chapter 2. The solution for this numerical problem is more complex than as presented here. Interested readers should consult ABAQUS theory manual [28]. 102 However, the use of CREEP user subroutine available in ABAQUS/Standard, as discussed in Section 4.2, has greatly simplified the problem to the end users whose only work is to define the relationship between hydrostatic and deviatoric stresses and the plastic strain rate. Adopting the Newtonian fluid model, our sintering model then has the form [54- 55]: 101.. 8'7 1 _ _ .. 81) -2# +3Kp( P 0051) (4.12) P All the symbols used in equation (4.12) are consistent with their notations as they appear in this thesis. It is important to note that when p = -as, and Sij = 0, then 35’ = O, as discussed in the Section 4.12. A look at equation (4.12) reveals a very important assumption of the sintering model: sintering is divided into two distinguishable modes, namely shear deformation and densification. The first part at the right of equation (4.12) deals with the earlier mode generally known as creep, whereas the second part the latter. Rahaman et al. [47] and Venkatachari and Raj [48] have verified experimentally the independence of the two modes. Furthermore, Rahaman et al. [47] verify experimentally that the creep rate increases linearly with applied shear stress. Venkatachari and Raj [48] on the other hand find that the densification rate is linearly proportional to the sum of the applied pressure and the sintering stress. The remaining question is how the shear and bulk Viscosities depend on 103 the most important parameters of sintering, such as porosity (or relative density), grain size, sintering temperature and sintering time, as outlined by Barsown [11]. 4.1.4 SHEAR AND BULK VISCOSITIES Bordia and Raj [79] presented a simple sintering model of ceramic film constrained by a rigid surface. This model follows the traditional assumption that sintering behavior can be separated into densification and shear modes, which are related directly to bulk and shear Viscosities respectively. The major distinction in this model is that the sintering behavior of the constrained ceramic film is dependent on the relative rates of shear and densification. Specifically, by assuming a simplified one-dimensional stress field and the applicability of phenomenological descriptions of the densification and shear properties, Bordia and Raj [79] obtained a time dependent analytical viscoelastic solution for the sintering behavior of the constrained ceramic film. The analytical solution shows that the shear stresses in a ceramic film may act as internal stresses that inhibit sintering if not relaxed by shear strains. The unrelaxed shear stresses may also result in tensile stresses within the film that may lead to sintering defects. This is an important insight of sintering behavior of films. However, there are still significant differences between the idealized model and our real compact. These include our real specimens being a compact of finite size and not a semifinite film, the varying material properties of our real specimens, and the nonlinearity in 104 the actual stress fields. While acknowledging their model, a different approach has been used in the current work. The following observations by Rahaman et al. [47] and Venkatachari and Raj [48] are definitive 'in determining the parameters dependence of shear and bulk Viscosities of the ceramic compact in this work: Venkatachari and Raj [48] propose for the densification rate to vary proportionally to the first power of the sintering stress, and creep rate to the first power of the effective stress. They also suggest that both densification and shear rates to be inversely proportional to the third power of the grain size [48]. Rahaman et al. [47] agree with their analysis of densification rates. With regard to shear, however, they maintain an exponential dependency of creep rate on density during the intermediate sintering stage when grain boundary diffusion is the operative mechanism. They believe that Coble’s relation of inverse dependency on relative density (3;; oc 71;) hold during the final sintering stage. This dependency has been traditionally described in terms of the stress concentration factor resulted from the ratio of the cross-sectional area of the specimen to the actual internal load bearing area. Rate of densification slows with increasing density due to decreasing sintering potential [47, 48]. The decrease in sintering pressure with increasing density may appear counterintuitive. As the pores become smaller, their curvature, and therefore the sintering pressure should 105 increase. However, in practice, this increase is often overcome by the rapid decrease in the number of pores. . Strain rate increases with increasing porosity but is sensitive to grain size. At a given grain size, however, a lower green density might experience lower strain rate because of higher grain growth [47]. . Applied load has no effect on relative density (or porosity) if the load is smaller than the sintering potential. Otherwise, it will increase the densification rate [48]. . Densification rate is almost independent of the initial compact density [47]. These observations are mostly met in Kwon et al.’s [55-57] two-stage sintering model: . 1- D 5 52 3C 1 D(D- Do) kTR ZDS 2(D- IDO) (4.13) during the initial and intermediate stage (D<0.85) in which _ 2 3 Kp = D(D DO)2kTR (4.14) 3C2(1-Do) 5.11119 D2(D—D0)kTR3 and )up = (4.15) 3Cl (1 - D0 )gefirz And 1 $2 pl _lgejf 3C3s Sij-‘I'C4a D)2(—p- 0'55) ij (4.16) 106 during the final stage (0.95 \-\_1 Hi L ~7KK~H % 8 1 ,5 8 ‘1 Hyperbolic Model Modified Drucker-Prager/cap Model Figure 5.10. Relative density distribution of compacted and finally sintered disc. 124 The final average relative density and final grain size can be computed using the following equations: Dave _ ZViéi _ZV ~ (5.3) (5.4) Table 5.1 shows the final average density and grain size for cases shown in Figure 5.6 and Figure 5.7 using base states from both compaction models. Apparently, the final average relative densities are very similar in all four cases, regardless of the geometry or the compaction model. Since the grain size model is affected only by sintering time and sintering temperature, which are identical for the cases with and without embedded graphite pieces, the resulting grain size should be equal in all the cases, as evident from the simulations. Table 5.1. Average relative density and grain sizes for different cases Case Model Relative density Grain size ( u m) Figure 5.6 hyperbolic 0.9608 0.4324 D-P/cap 0.9599 0.4324 Figure 5.7 hyperbolic 0.9568 0.4324 D-P/cap 0.9552 0.4324 As a final note, the effect of elastic strain is assumed negligible in the densification rate equation (4.20). This assumption is verified by comparing the 125 compact mass before and after final sintering. For all the simulated cases, the mass difference before and after sintering simulations is at most 0.2%, justifying the negligible elastic strain effects. Thus, the densification rate equation is consistent with the principle of mass conservation. 5.3.2 SINTERING TIME One of the goals in sintering is to achieve a relative density close to the theoretical values at a reasonably low temperature in the shortest time. If an accurate material data is available, the prediction of sintering time (only the isothermal heating duration, not the heating up and cooling down periods) to achieve desired relative density is highly probable. This will reduce manufacturing cost and time. The specimens in our lab do not sinter to theoretical density at 1475°C for 4 hours in a conventional furnace. However, zirconia tends to sinter up to a relative density of 0.95 in our experiment. It is observed from the sintering simulations that the effects of compaction model and geometry on sintering time are not significant. Since the final relative density distribution has been shown to be quite uniform, the sintering time needed to achieve a certain relative density in a compact with complex geometry may be obtained from that of a simple geometry. 126 5.4 CONCLUSIONS The sintering model has shown its capability in capturing the final sintered shape of the compacts, the final relative density and residual stresses, and the required sintering time to achieve a desired densification. The model predicts uniform relative density distribution and grain size regardless of the compaction model or compact geometry. To refine this model, a sinter forging experiment has to be performed to measure experimental parameters C1 and Ca. The fine-tuned model will then be applied to simulate the micro heat exchanger fabrication process. 127 Chapter 6 CONCLUSIONS 6.1 CONCLUSIONS A finite element model for compaction and sintering of ceramic powder based on plasticity theory are presented. The compaction model has the capability to predict the relative density distribution and the residual stress in a green compact. It also predicts the location of microcracks based upon regions of shear yielding during the top punch removal step. The sintering model predicts the shrinkage of the green compact and grain growth during sintering. Material data for the current compaction simulations have been adopted from open literature and experiments. The simulation results are compared to the experimental results available in the literature and found to be in a reasonable agreement. Most material parameters for the current sintering simulations have been curve- fitted from data obtained from pressureless sintering experiments. The goal of the current research is to apply these models to the micro heat exchanger in development. Though current experimental findings are believed to be consistent with the simulation results, further efforts are needed to improve the accuracy of the model. These efforts are discussed in the following section, and most of them require facilities currently not available in our institution. 128 0" 6.2 FUTURE WORKS The implicit assumption in comparing our simulation results with our own experimental findings is that the powder used possesses a similar mechanical behavior to those in the literature. This assumption needs to be verified. We may either perform tests on the powder we use, compare our material data to those in the literature and then redo our simulations with the new material parameters; or use the same powder whose material parameters were determined in the literature. Our compaction simulation results indicate that the modified Drucker-Prager/cap model and the hyperbolic model differ in predictions only when the applied axial pressure exceeds a certain limit. The simulation results also suggest that the differences become increasingly obvious beginning from the top punch removal step. This leads us to believe that the unloading of the powder compact plays a more significant role than expected. An in-depth investigation needs to be conducted to better model the compact behavior during stress relaxation. This will inevitably require shear yield testing on the powder, and an experimental methodology to quantify the yielding process. Development along this line may help establish a critical shear yield strain value upon which micro-cracks will occur. The limitation of our lab facilities has deterred efforts in this direction. 129 To better model the sintering behavior, experimentation on sinter forging needs to be conducted to quantify shear strain rates of the powders used. These experiments will yield the data for the improvement in predicting the shape distortions of the sintered products. Sinter forging has to be conducted with at least a dilatometer, a piece of equipment we currently do not possess. Having partial material data from pressureless sintering and not sinter forging may be insufficient to justify the current model. The experimental data collected in the suggested compaction and sintering tests will improve the quality of the material data needed in our simulations. However, it is uncertain whether this improvement will be drastic or even feasible. On the other hand, a more direct correlation of experimental findings with the simulation results may be difficult to establish without the suggested improvement. Finally, the current simulations are performed only on zirconia data. Sintering data on alumina is available in literature, but the compaction data are not available in the low relative density range. The relative density of the alumina powder required in the experiment is only about 0.23, and most of the ones available in the literature are in the range of 0.6. Extrapolations on the data are attempted but they have not been fruitful and lack justifications. 130 APPENDIX A: Partial Derivatives for Hyperbolic Cap Surface The following are the partial derivatives needed for Newton’s iterative method in integrating the elastic-plastic equations with hyperbolic yield function f1(p, q,D)= O=—+cosh{A(D)l p —}- B(D). 106106 af 1 _=_ A.1 34 106 ( ) 3f: A—I—lg)sirm h{A(D) ”—6} (A2) p 10 106 _f_:0 3p aq (A.3) _&)f__ p £171__(_D)$i p 33(0) 30 10—6 30 nh{A(D)10—6} 30 (AA) 32f _ 1 34(0)si _p_ p Ame/11mm p 31230 106 30 nh{A(D)__106 106 106 3D hA(D)1—06 (A.5) 32f _ 37133-0 (A.6) 1M) 1 A) _= __ A(D) A.7 ap2 106 106 ( ) 2 §_f=o (A.8) aqz 131 APPENDIX 8: Newton Iteration for Compaction Constitutive Model The right-hand—side of equations (2.38) and (2.39) can be expanded by partial differentiation such that 2 2 82 13135 +1»: a—{awa faq+a 82%)] +afdA£q +135 [apfiaw Hag—2130 =0 34 p p 31734 aq2 813961 3P <3an 3—WD and (8.1) flaf-aeralanr-aiao =0 (3.2) 3p dq GD where 3p = ape +d(KA£p) = —K[ : a§+ KdAep aq = dqe + 3(-3/1A£ p) = 2,113 : dg- 3M3A£p (3.3) 8h 1 ah ah ah ah GD: l—— 8A ——8A —8 —3 [ an) [amp 8" 321 Age 8" Nap “aqq ‘1] In (3.3), the moduli K and G are assumed to be constants. During each iteration to solve for Asp and Aeq, the strain states (and the relative density at the beginning of the increment) are assumed constant. Therefore, the solution method can be thought as incremental linear approach, where g can be assumed constant. Apparently, in solving for the linearization moduli, the same assumption cannot be applied since g vary from increment to increment. 132 Implementing equations (2.40), (2.41) and noting equations (B.1)-(B.3), coefficients A11,A12,A21,A22,b1, b2 can be shown to be af 82f 82f an 32f 32f an A =— A K A K ‘1 aq+ 8 6[ apaq+dan 6A8 + £6 apz +apap BASP af 82f 82f an 32f 82f an A =— A 3 A5 -3 ‘2 ap+ £P[ ”3,123an 3A8 + 6 ”am 3,030 3qu 91 a: 80 A21: ap an amp (3'4) af 3f 30 A =- — — 3”aq aDaAeq Expressions for partial differentials can be found in Appendix A for yield function (2.14a). Explanation on (8.30) Adopting Aravas [72] to fit our current one parameter model, an ( ah j" ah ah = 1- — + K— BAEP BD BAEP 8p (3.5) an =[1— a_h_) [ah _3#8_h] dAe‘q 8D dAeq aq In our case, 1') = —Dé‘kk =h{ép ,é‘q,p, q,D} (3.6) AD = —DA£kk= {ASP ,Ab‘q , p,q,D} From (8.6), a—D = 4,5: 133 But as observed from (2.30), Asp = A85, ‘ .p _ . In other words, Ekk — 3,, We can now write (8.6) as AD = —DAE,5c = -D(A€p )= h{A£p ,Aeq,p,q,D} It is now obvious from (2.44) that ah = —D, % = O, 5% = 0 dAep ap aq Equations (3.5) then reduces to DD aAsp 3D aAsq = + 315‘ )—1(- D) 134 (13.7) (8.8) (3.9) (8.10) (3.11) fl'r APPENDIX C1: Linearization Moduli for Shear Yield Surface In the following development, all the symbols, superscripts and subscripts used are consistent with definitions used in this work thus far. The deviatoric strain is g=§~1trace(§>z=[1—lu)§=g6 +g” (01.1) 3 3 T- Recall that §t+At = zflgfflv = 2fl(§te + A§t+At - Aim) (C12) Now we define the equivalent deviatoric strain increment A?" such that Agp = AEPQ (01.3) Combining (C12) and (CI .3) yields (IMA: where §=gf +A§1+At (C1.5) Along the critical state line, non-dilational deformation with the proposed non- associative flow rule, as summarized in equations (2.27) and (2.28), dictates the plastic strain rate to be independent of hydrostatic stress, p. In other words, the variation of (01.4) with respect to p is zero. 135 Now, taking the variation of (01.4) with respect to all the other quantities at the end of the increment gives 3 -p 3 -p A2” ,. [l+-—’u—Aet+m]d§,+m+§t+m ,u [3A6 - dqt+m]=2,uag 41+At qr+At qt+At (01.6) But taking the inner product of (01.4) to itself, q" = cm ”#125,, =3”? (01.7) .. 2 . . where e = Egzg (01.8) Again, taking the variation of (01.7) yields (3qu + Max:211,“ =3uaz; (01.9) Assumption of zero hardening along this line translates to zero variation of q with respect to equivalent plastic strain. Therefore, (C1.9) reduces to 3112:,“ = 33%, (01.10) From the definition in (01.8), az=3zgzag ((21.11) 3e Combining (C1.10) and (01.11) -p .. 2 .. .. aAeHAt =86 zg-Egzag (C1.12) Substituting (01.12) into (C1.6), and after some manipulation, a§t+At :[Q"1_R”§t+m§.t+m]:aé (C113) where 136 Q"-_-_2.qt:At =2” qt+Af (C1.14) 3 e qe R": l — 3” (01.15) qt+AtE Qt+Ath Since it is obvious from (017) that 2: = 512 (01.16) 3;! Recognizing that at+At :St+At "PHAtl (C1-17) aé=(_J_-%I_I.]3§,+At (01.18) and because of plastic flow being independent of hydrostatic pressure, p 1 BpHA, = -KL : a; (C1.19) finally, 91"” =Q"1+ (K —-;—Q")u-R"g,+m §,+A, (01.20) 137 APPENDIX C2: Linearization Moduli for Hyperbolic Cap Surface Since g = Zp(—%tr(§); + g), (02.1) the deviatoric trial stress can be written as g, = s + 241 - éfljAg (02.2) Consequently, age = 22(1 - 111) (02.3) a§t+At 3 Rewriting equation (2.6), we have qe2 =2s 1s (02.4) 2 — e -e Differentiating (02.4), aqe _ aqe age _3: a§e — ‘ = 2 02.5 3; age 3; 24c 32%, Mum ( ) - 3 8m” ””6” =2_-‘-e and 6m!” (02.6) e So 621m: 3 32., _ 3 acre =2_#[§J_1”_n n ] a§t+At zqe a£t+m 2qe2 39M, qe 2“ 2— -H'A! —t+At (02.7) Using equations (8.1), (8.2) and (8.3) in Appendix 8, we arrived at: C11=A21. C12=A22. 021=A11. 022=A12. (CZ-8) where the A); coefficients are those in Appendix 8, 138 and 12.2 a_f( a2) a2) ,0... 6p 8D 8D 8p f 3f ah J I ah D - -—2 — C2.9b 12 ”[34. aq 8D __l—( 3D dq ( ) 32f 32f a2 f a_2f ( 8h)_13h D =K A — A —— K A A l-— 2‘ [ £6 apaq+6 86 3,,2 2)+ [ £6 aoaq+ 86— map an ap (02.9c) r 82f 32f 32f 32f 91th D = -2 A — A 2 A A — l- — 22 ,u[ 8” Maq + 8" aqdp] G[ 8” £9qu + sq aDap 3D aq (02.9d) 139 APPENDIX D: Constitutive Laws for lncompressible Fluid Constitutive laws for incompressible fluid can be summarized in the equations below. 61' = 266:7 ”(9:32, . . l 1 81=£f +317—[01 -§(0'2 +03):l . . 1 l 82 = 8f +-3;I:0'2 -§(O'1+O'3):| . . l 1 £3 = sf +—3—5[a3 -§(01+02)] 140 (0.1) (0.2) (D.3a) (D.3b) (D.3c) BIBLIOGRAPHY 1 Drucker, D. 0., and W. Prager, “Soil Mechanics and Plastic Analysis or Limit Design,” Quarterly of Applied Mathematics, Vol. 10, 1952, pp. 157-165. 2 Chen, W. F. and D. J. Han, Plasticiy for Structural Engineers, 1988, New York: Springer-Vedag. 3 Schofield, A., and P. Wroth, Critical State Soil Mechanics, 1968, London: McGraw Hill. 4 Roscoe, K. H., A. N. Schofield, and A. Thurairajah, “Yielding of Clays in States Wetter than Critical,” Geotechnique, Vol. 13, 1965, pp. 211-240. 5 Hibbitt, Karlsson and Sorensen, Abagus User Manual Volme l Vprsion 5.7, 1998, Rhode Island, Hibbitt, Karlsson and Sorensen, Inc. 6 Bortzmeyer, D., “Modeling Ceramic Powder Compaction,” Powder Technology, Vol 70, 1993, pp. 131-139. 7 Aydin, I., 8. J. Briscoe and K.Y. Sanliturk, “The lntemal Form of Compacted Ceramic Components: A Comparison of A Finite Element Modeling With Experiment,” Powder Technology, Vol. 89, 1996, pp. 239-254. 8 Aydin, I., 8. J. Briscoe and K. Y. Sanliturk, “Dimensional Variation of Die- Pressed Ceramic Green Compacts: Comparison of a Finite Element Modelling with Experiment,” Journal of European Ceramic Society, Vol. 17, 1997, pp. 1201- 1212. 9 Kim, K. T., S. W. Choi and H. Park, “Densification Behavior of Ceramic Powder Under Cold Compaction,” Journal of Engineering Materials and Technology, Vol. 122, 2000. PP. 238-244. 10 Zipse, H., “Finite-Element Simulation of the Die Pressing and Sintering of a Ceramic Component,” Journal of European Ceramic Society, Vol. 17, 1997, pp. 1707-1713. 11 Barsown, M., “Chapter 10: Sintering and Grain Growth,” Fundamental of Ceramics, New York: The McGraw-Hill Companies, Inc., 1997, pp. 331 -385. 12 Arzt, E., “The Influence of An Increasing Particle Coordination on the Densification of Spherical Powder,” Acta Metallurgica, Vol. 30, 1982, pp. 1883- 1890. 141 13 Fischmeister, H. F., E. Arzt, and L. R. Olsson, “Particle Deformation and Sliding During Compaction of Spherical Powders: A Study by Quantitative Metallography,” Powder Metallurgy, Vol. 21, No. 4, 1978, pp. 179-186. 14 Fischmeister, H. F. and E. Arzt, ”Densifciation of Powders by Particle Deformation,” Powder Metallurgica, Vol. 26, 1983, pp. 82-88. 15 Artz, E., M. F. Ashby and K. E. Easterling, “Practical Application of Hot- lsostatic Pressing Diagrams: Four Case Studies,” Metallurgica Transactions A, Vol. 14A, February 1983, pp. 211-221. 16 Li, E. K. H., and P. D. Funkenbusch, “Experimental Size Ratio and Compositional Effects on the Packing and Hot lsostatic Pressing of Spherical Powders,” Materials Science and Engineering, A157, 1992, pp. 217-224. 17 Li, E. K. H. and P. D. Funkenbusch, “Modeling of the Densification Rates of Monosized and 8imodaI-Sized Particle Systems During Hot lsostatic Pressing (HlP),” Acta Metallurgica, Vol. 37, No. 6, pp. 1645-1655. 18 Nair, S. V., 8. C. Hendrix, and J. K. Tien, “Obtaining the Radial Distribution of Random Dense Packed Hard Spheres,” Acta Metallurgica, 1986, pp. 1599-1605. 19 Nair, S. V. and J. K. Tien, “Densification Mechanism Maps for Hot lsostatic Pressing (HIP) of Unequal Sized Particles,” Metallurgical Transactions A, Vol. 18A, January 1987, pp. 97-107. 20 Aydin, |., 8. J. Briscoe and N. Ozkan, “Modeling of Powder Compaction: A Review,” MRS Bulletin, Vol. 22, n. 12, DEC 1997, pp. 45-51. 21 Chen, W. F., “Chapter 5: Constitutive Relations for Concrete, Rock and Soils: Discusser’s Report, Section II,” Mechanics of Geomaterials: Rocks, Concrete, Soils, ed. 2. P. Bazant, 1984, Chichester: John Wiley & Sons Ltd, pp. 74-86. 22 Walker, D. M., “An Approximate Theory for Pressures and Arching in Hoppers,” Chemical Engineering Science, Vol. 21, 1966, pp. 975-997. 23 Briscoe, 8. J. and N. Ozkan, “Compaction Behavior of Agglomerated Alumina Powders,” Powder Technology, Vol. 90, 1997, pp. 195-203. 24 Chen, W. F., “Chapter 5: Constitutive Modeling in Soil Mechanics,” Mechanics of Engineering Materials, ed. Desai, 0. S, and Gallagher, R. H, 1984, Chichester: John Wiley & Sons Ltd, pp. 91-120. 25 Aydin, l., 8. J. Briscoe, and K. Y. Sanliturk, “Density Distributions during the Compaction of Alumina Powders: A Comparison of A Computation Prediction With Experiment,” Computational Materials Science, Vol. 3, 1994, pp. 55-68. 142 26 Simo, J. C. and T. J. R. Hughes, “Computational lnelasticity,” lnterdisplinary Applied Mathematics, Vol. 7, 1998, New York: Springer. 27 Shima, S. and M. Oyane, “Plasticity Theory For Porous Metals,” lntemational Journal of Mechanical Science, Vol. 18, 1976, pp. 285-291. 28 Hibbitt, Karlsson and Sorensen, Abagus Theorv Manual Version 5.8, 1998, Rhode Island: Hibbitt, Karlsson and Sorensen, Inc. 29 Chen, W. F., “Soil Mechanics, Plasticity and Landslides,” Mechanics of Material Behavior: The Drucker-Prager Anniversary Volume, ed. George J. Dvorak and Richard T. Shield, 1984, Amsterdam: Elsevier, pp. 31-58. 30 Mori K., ”Finite Element Simulation of Nonuniform Shrinkage in Sintering of Ceramic Powder Compact,” ln: Chenot JL et al., editors, Mmerical Methods in Industrial Forming Processes NUMIFORM 92, Rotterdam: Balkema, 1992. pp. 69-78. 31 Mori, K-I, K. Osakada and T. Hirano, “Finite Element Modelling of Nonuniform Shrinkage In Sintering of Ceramic Powder Compact,” Computer Aided Innovation of New Materials II, 1993, pp. 1781. 32 Mori, K-I., K. Osakada , and M. Miyazaki, “Prediction of Fracture in Sintering of Ceramic Powder Compact,” International Journal of Machine Tools Manufacture, Vol. 37, No. 9, 1997, pp. 1327-1336. 33 Mori K, Osakada K., “Net Shape Approach for Sintering Process of Graded Laminated Powder Materials using Finite Element Simulation,” Annals of the CIRP Vol. 48, No. 1, 1999, pp. 239-242. 34 Shima, S. and K. Mimura, “Densification Behavior of Ceramic Powder,” lntemational Journal of Mechanical Sciences, Vol. 28, No. 1, 1986, pp. 53-59. 35 Yeh, T. S., and M. D. Sacks, “Effect of Particle Size Distribution on Sintering of Alumina,” Journal of American Ceramic Society, Vol. 71, No. 12, 1988, pp. 0- 484-0-487. 36 Ting, J. M. and R. Y. Lin, “Effect of Particle-Size Distribution On Sintering Part I: Modelling,” Journal of Materials Science, Vol. 29, 1994, pp. 1867-1872. 37 Ting, J. M. and R. Y. Lin, “Effect of Particle-Size Distribution On Sintering Part II: Sintering of Alumina,” Journal of Materials Science, Vol. 30, No. 9, 1995, pp. 1867-1872. 143 38 Darcovich, K., L. Bera, K. Shinagawa, “Particle Size Distribution Effects in an FEM Model of Sintering Porous Ceramics,” Material Science and Engineering, A341, 2003, pp. 247-255. 39 Darcovich, K., T. Floyd, P. Hontanx, V. Roux and K. Shinagawa, “An Experimental and Numerical Study of Particle Siaze Distribution Effects on Sintering of Porous Ceramics,” Material Science and Engineering, A348, 2003, pp.76-83. 40 Swinkels, F. B. and M. F. Ashby, “Overview 11: A Second Report On Sintering Diagrams,” Acta Metallurgica, Vol. 29, 1980, pp. 259-281. 41 Lippmann, H., and R. Iankov, “Mathematical Modeling of Sintering During Powder Forming Processes,” lntemational Journal of Mechanical Science, Vol. 39, No.5 1997, pp. 585-596. 42 Bordia, R. K., and G. W. Scherer, “On Constrained Sintering—l. Constitutive Model For A Sintering Body,” Acta Metallugica, Vol. 36, No. 9, 1988, pp. 2393- 2397. 43 80rdia, R. K. and G. W. Scherer, “On Constrained Sintering—II. Comparison of Constitutive Model,” Acta Metallugica, Vol. 36, No. 9, 1988, pp. 2399-2409. 44 Raj, R. and R. K. Bordia, “Sintering Behavior of Bi-modal Powder Compacts,” Acta Metallurgica, Vol. 32, No.7, 1984. PP. 1003-1019. 45 Hsueh, C. H., A. G. Evans, and R. M. Cannon, “Viscoelastic Stresses and Sintering Damage In Heterogeneous Powder Compacts,” Acta Metallurgica, Vol. 34, No. 5, 1986, pp. 927-936. 46 Scherer, G. W., “Sintering Inhomogeneous Glasses: Application To Optical Waveguides,” Journal of Non-Crystalline Solids, Vol. 34, 1979, pp.239-256. 47 Rahaman, M. N., L. C. De Jonghe, and R. J. Brook, “Effect of Shear Stress on Sintering,” Journal of American Ceramic Society, Vol. 69, No. 1, 1986, pp. 53-58. 48 Venkatachari, K., and R. Raj, “Shear Deformation and Densification of Powder Compacts,” Journal of American Ceramic Society, Vol. 69, No. 6, 1986, pp. 499-506. 49 Riedel, H. and D.-Z. Sun, “Simulation of Die Pressing and Sintering of Powder Metals, Hard Metals and Ceramics," Numerical Methods in Industrial Forming Processes, ed. Chenot, Wood and Zienkiewicz, 1992, Rotterdam: Balkema, pp. 883-886. 144 50 Du, Z.-Z. and A. C. F. Cocks, “Constitutive Models for the Sintering of Ceramic Components—I. Material Models,” Acta Metallurgica, Vol. 40, No. 8, 1992, pp. 1969-1979. 51 Du, Z.-Z. and A. C. F. Cocks, “Constitutive Models for the Sintering of Ceramic Components—II. Sintering if Inhomogeneous Bodies,” Acta Metallurgica, Vol. 40, No. 8, 1992, pp. 1981-1994. 52 Shinagawa, K., “Micromechanical Modelling of Viscous Sintering and A Constitutive Equation with Sintering Stress,” Computational Materials Science, Vol. 13, 1999, pp. 276-285. 53 Kraft, T. and H. Riedel, “Numerical Simulation of Die Compaction and Sintering,” Powder Metallurgy. Vol. 45, No. 3, 2002, pp. 227-231. 54 Kraft, T., H. Riedel, P. Stingl, and F. Wittig, “Finite Element Simulation of Die Pressing and Sintering,” Advanced Engineering Materials, Vol. 1, No. 2, 1999, pp. 107-109. 55 Kwon, Y. S., G. Son, J. Suh and K. T. Kim, “Densification and Grain Growth of Porous Alumina Compacts,” Journal of American Ceramic Society, Vol. 77, No. 12, 1994. PP. 3137-3141. 56 Kwon, Y. S. and K. T. Kim, “High Temperature Densification Forming of Alumina Powder—Constitutive Model and Experiments,” Journal of Engineering Materials and Technology, Vol. 118, 1996, pp. 448-455. 57 Kwon, Y. S. and K. T. Kim, “Densification Forming of Alumina Powder— Effects of Powder Law Creep and Friction,” Journal of Engineering Materials and Technology, Vol. 118, 1996, pp. 471-477. 58 Kim, K. T., Y. S. Kwon, and H. G. Kim, “Near-Net-Shape Forming of Alumina Powder under Hot Pressing and Hot lsostatic Pressing,” lntemational Journal of Mechanical Science, Vol. 39, No. 9, 1997, pp. 1011-1022. 59 Kim, K. T, H. G. Kim, and H. M. Jang, “Densification Behavior and Grain Growth of Zirconia Powder Compact under High Temperature,” lntemational Journal of Engineering Science, Vol.36, 1998, pp. 1295-1312. 60 Kim, H. G., H. M. Lee, and K. T. Kim, “Near-Net-Shape Forming of Ceramic Powder under Cold Combination Pressing and Pressureless Sintering,” Journal of Engineering Materials and Technology, Vol. 123, 2001, pp. 221-228. 61 Cocks, A. C. F., “The Structure of Constitutive Laws for the Sintering of Fine Grained Materials, ” Acta Metallurgica, Vol. 42, No. 7, 1994, pp. 2191-2210. 145 62 Kim, H. G., O. Gillia and D. Bouvard, “A Phenomenological Constitutive Model for the Sintering of Alumina Powder,” Journal of the European Ceramic Society, Vol. 23, 2003, pp. 1675-1685. 63 Kim, H. G., O. Gillia, P. Doremus and D. Bouvard, “Near Net Shape Processing of A Sintered Alumina Component: Adjustment of Pressing Parameters Through Finite Element Simulation,” lntemational Journal Mechanical Sciences, Vol. 44, 2002, pp. 2523-2539. 64 Jagota, A., K. R. Miskeska, and R. K. Bordia, “Isotropic Constitutive Model for Sintering Particle Packings,” Journal of American Ceramic Societym Vol. 73, No. 8, 1990. PP. 2266-2273. I 65 Sanliturk, K. Y., I. Aydin, and 8. J. Briscoe, “ A finite-Element Approach for the Shape Predition of Ceramic Compacts during Sintering,” Journal of American Ceramic Society, Vol. 82, No. 7, 1999, pp. 1748-1756. 66 Raj, R. and M. F. Ashby, “Intergranular Fracture at Elevated Temperature,” Acta Metallurgica, Vol.23, 1975, pp. 653-666. 67 Bouvard, D., P. Doremus, O. Gillia and V. Bomnefoy, “Finite Element Simulation of Compaction and Sintering of Ceramic Powders,” Key Engineering Materials, Vol. 206-213, 2002, pp. 243-248. 68 Zeuch, D. H, J. M. Grazier, J. G. Arguello, K. G. Ewsuk, “Mechanical Properties and Shear Failure Surfaces for Two Alumina Powders in Triaxial Compression,”Journal of Materials Science, Vol. 36, 2001, pp. 2911-2924. 69 Westerheid, R., Zipse, H. and Hollstein, T. “Finite Element Simulation of Die Pressing and Sintering Results in Time and Cost Savings,”0eramic Forum lntemational (CFl)/8erichte der DKG, Vol. 75, No. 4, 1998, pp. 31 -34. 70 Briscoe, 8. J. and Rough, S. L., “The Effects of Wall Friction in Powder Compaction,” Colloids and Surfaces A: Physicochemical and Engineering Aspects, Vol. 137, 1998, pp. 103-116. 71 Govindarajan, R. M., Deformation Processing of Porous Metals, Doctoral Thesis, 1992, University of Pennsylvania. 72 Aravas, N., “On the Numerical Integration of a Class of Pressure-Dependent Plasticity Models,” lntemational Journal of Numerical Methods in Engineering, Vol. 24, 1987, pp. 1395-1416. 73 Loret, 8. and Prevost, J. H., “Accurate Numerical Solutions for Drucker-Prager Elastic-Plastic Models,” Computer Methods in Applied Mechanics and Engineering, Vol. 54, 1986, pp. 259-277. 146 74 Bonet, Javier and Richard D. Wood, Nonlinear Contimym Mecha_nics For Finite Element Analysis, 1997, Melbourne: Cambridge Univerity Press. 75 Shin, H., Kok, 0., Case, 8, Kwon, P., (2003), “Meso—Scale Channels On Ceramics 8y Machining,” ASPE Winter Topical Meeting, Machines and Processes for Micro-Scale And Mesa-Scale Fabrication, Metrology and Assembly, Gainesville, Florida, in publication. 76 Shin, H., Case, 8, Kwon, P., (2002), “Novel Powder Processing Techniques to Fabricate an Efficient Meso-scale Heat Exchanger,” Transactions of NAMFIII, XXX. PP. 671-679. 77 Helle, A. S., K. E. Easterling and M. F. Ashby, “Hot-isostatic Pressing Diagrams: New Developments,” Acta Metallurgica, Vol. 33, No. 12, 1983, pp. 2163-2174. 78 McMeeking, R. M. and L. T. Kuhn, “A Diffusional Creep Law for Powder Compacts,” Acta Metallurgica, Vol. 40, No. 5, 1992, pp. 961-969. 79 Bordia, R. K. and R. Raj, “Sintering Behavior of Ceramic Films Constrained by a Rigid Substrate,” Journal of American Ceramic Society, Vol. 68, No. 6, 1985, pp. 287- 292. 147 MMMMMMM Milli/ill llillllllilllllllellil 93 02504 26 4 ‘2 ~ g...—-—._. _ - v f...— _—