Thermomechanical modelling of ceramic pressing and subsequent sintering
TThermomechanical modelling of ceramic pressing andsubsequent sintering
D. Kempen, A. Piccolroaz, and D. Bigoni ∗ Department of Civil, Environmental and Mechanical Engineering, University of Trento, Italy
Abstract
An elastic-visco-plastic thermomechanical model for the simulation of cold form-ing and subsequent sintering of ceramic powders is introduced and based on mi-cromechanical modelling of the compaction process of granulates. Micromechanicsleads to an upper-bound estimate of the compaction curve of a granular material,which compares well with other models and finite element simulations. The param-eters of the thermomechanical model are determined on the basis of available dataand dilatometer experiments. Finally, after computer implementation, validation ofthe model is performed against experiments developed on specially designed ceramicpieces, characterized by zones of different density. The mechanical model is found toaccurately describe forming and sintering of stoneware ceramics and can therefore beused to analyze and optimize industrial processes involving compaction of powdersand subsequent firing of the greens.
The production of ceramic pieces is based on technologies involving a massive waste ofenergy and materials , so that environmental preservation imposes a rationalization ofindustrial processes to reduce pollution. The optimization of the production process isdirectly linked to the availability of models for the mechanical behaviour of the powdersand binders used during forming, including the simulation of both phases of cold com-paction and subsequent sintering, and in the design of mechanical characteristics of thefinal pieces.Sintering is the common method for completing the process of ceramic production andinvolves heating of the green bodies, obtained in a preceding compaction stage, to obtainthe required density and strength of the final products. Understanding and modelling themechanisms involved in the processes of compaction and sintering is therefore crucial toensure high quality and reproducibility of ceramic materials. ∗ Corresponding author: e-mail: [email protected]; phone: +39 0461 282507. Grinding of the raw material requires mills with a power up to 1 MW, while drying and increasingof the temperature of the slip involves up to 500 KW of electrical power and 15 Gcal/h of thermal power.The 80 % of the thermal capacity is lost at the chimney and powder is spread in the environment. Theforming of the ceramic powders represents a huge waste of energy, because only the 5 % of the energy istransmitted to the final piece from the presses (up to 250 KW of installed power). Finally, drying andsintering requires large burners consuming up to 10 Gcal/h. a r X i v : . [ c s . C E ] A p r n the last 50 years, there have been significant developments in the theory of sinter-ing, but only in the 1980s the use of continuum mechanics was introduced, to predict thestresses and strains that can develop during the process. A comprehensive review of themodels for sintering proposed in the last 25 years can be found in [1]. In recent years,meso scale models have been developed using multi particle Finite Element (MPFEM)approaches [2] or discrete element approaches [3]. These models are useful in the under-standing of micro- or meso-mechanical effects, but, due to the limitation in the numberof particles that can analyzed, they are hardly applicable to model whole complex shapedparts, so that continuum models, such those developed in [4–6], play a central role.Often sintering was studied with reference to isothermal conditions [6, 7], but duringthe firing the body is subject to non-isothermal loadings, especially for high heating ratesand large ceramic pieces, so that the modelling of this situation is one of the objectivesof the present article.Before sintering, a green body is obtained through powder compaction and, depend-ing on the geometry and the compaction method, the pressed green body shows usuallysignificant differences of local density throughout the part [8, 9], so that powder pressinghas a strong influence on sintering, which depends on the (relative) density of the green.Therefore a continuum model that can predict the density distribution after the pressingstep and the subsequent sintering is important for industrial applications [10]. In thisresearch direction, while cold compaction has been often addressed [11–19], the combina-tion of pressing and sintering has been scarcely investigated. In particular, two differentmodels, one for powder pressing and another for sintering have been proposed [7], whilethermomechanical models for hot-isostatic pressing of metal powders were developed usinga pressure-sensitive and temperature dependent yield function [20, 21].The objective of the present article is to propose a thermomechanical, elastic-visco-plastic model to simulate cold powder compaction and subsequent non-isothermal solidstate sintering. The model is grounded on the thermomechanical framework developedby Simo and Miehe [22, 23] and on the additive decomposition of elastic and plastic log-arithmic strains [23, 24]. Moreover, the BP yield function [25] is used to simulate bothcompaction and firing, so that these two process steps become more closely integrated.With reference to an aluminum silicate spray dried powder (used for the industrial pro-duction of ceramic tiles), the model is calibrated against both available data and ad-hoc performed experiments and is implemented (through an UMAT routine) in a finite ele-ment code (Abaqus). Finally, validation against experiments, conducted on large ceramicpieces (in which density variations were intentionally introduced), shows that the modelprovides an accurate prediction of the entire process of forming and sintering and thereforeremains now available for the design and optimization of ceramic pieces. With the purpose of introducing a model for sintering of green ceramics, the large defor-mation of a thermo-elastic-visco-plastic solid is considered, described by the deformationgradient F , the right Cauchy-Green tensor and the corresponding Lagrangean logarithmicstrain C = F T F , (cid:15) = 12 log C , (1)where the superscript T denotes the transpose.2 key assumption is the additive decomposition of logarithmic strain into an elastic(subscript ‘ e ’) and a visco-plastic (subscript ‘ p ’) component proposed by Miehe et al. [23]and Sansour and Wagner [24] as (cid:15) = (cid:15) e + (cid:15) p . (2)The Helmholtz free-energy ψ ( (cid:15) e , T, (cid:98) ρ ) is introduced as a function of (i.) the elasticpart of the logarithmic strain, (ii.) the temperature T , and (iii.) the inelastic relativedensity of the material (cid:98) ρ , which is the only internal variable introduced in the treatment,defined as (cid:98) ρ = ρ ρ fd e − tr (cid:15) p , (3)and representing a dimensionless measure of the mass density upon unloading. In Eq. (3) ρ is the initial value of the mass density ρ , and ρ fd the value corresponding to the fullydense material, which contains no pores.Definition (3) may be better appreciated by writing the rate of mass conservation ˙ ρρ = − tr ˙ (cid:15) , (4)which can be integrated in time to provide the expression log ρρ = − tr (cid:15) . (5)From eq. (5) the rate of mass density can be calculated to be ˙ ρρ fd = [ − (cid:98) ρ tr ˙ (cid:15) e + ( (cid:98) ρ ) · ] e − tr (cid:15) e , (6)where ( (cid:98) ρ ) · = − (cid:98) ρ tr ˙ (cid:15) p , (7)an equation which will become useful later.Note also that the mass density is related to the porosity f (the ratio between thevolume of the voids and the total volume of a sample) of the material through the equation ρρ fd = 1 − f. (8)Eq. (8) can be understood considering the deformation gradient F fd and its determinant J fd needed to bring the current volume element V to the volume of the fully dense material V fd . In this deformation, the volume and the density transform according to the well-known rules V fd = J fd V and ρ = J fd ρ fd , respectively, so that eq. (8) is obtained, becauseby definition f = 1 − V fd /V .The so-called elastoplastic coupling (in which the plastic strain influences the elasticstiffness, [26–28]) is not introduced in the model and the elastic part of the deformation(not particularly important during sintering) will be eventually treated with the standardlinear isotropic thermoelastic law [29]. Therefore, the stress σ , work-conjugate in the Hillsense to the Lagrangian logarithmic strain, the thermodynamical force R associated tothe internal variable (cid:98) ρ , and the entropy η can be expressed as σ = ρ ∂ψ ( (cid:15) e , T, (cid:98) ρ ) ∂ (cid:15) e , R = − ρ ∂ψ ( (cid:15) e , T, (cid:98) ρ ) ∂ (cid:98) ρ , η = − ∂ψ ( (cid:15) e , T, (cid:98) ρ ) ∂T . (9)3or simplicity, an additive form of the Helmholtz free energy is assumed, sum of theelastic and the purely thermal energies, plus a ‘pore energy’ ψ ( (cid:15) e , T, (cid:98) ρ ) = ψ e ( (cid:15) e , T ) + ψ T ( T ) + ψ pore ( (cid:98) ρ ) , (10)where ψ T ( T ) = − c h T log TT + c h ( T − T ) − η ( T − T ) + ψ , (11)in which T is the absolute temperature corresponding to the unstressed material when (cid:15) = 0 ; moreover, η and ψ are the values of entropy η and free energy at T = T and (cid:15) = 0 ; finally c h = − T ∂ ψ ( (cid:15) e , T, (cid:98) ρ ) ∂T , (12)is the specific heat at constant values of strain and internal variables.The density, together with the (visco-)plastic strain (cid:15) p are internal variables of thesystem. The external variables are the total strain (cid:15) and the temperature T . To followthe theory of thermodynamics of irreversible processes, the Clausius-Duhem dissipationinequality has to be fulfilled at all times [29, 30], namely − ρ (cid:16) ˙ ψ + η ˙ T (cid:17) + σ · ˙ (cid:15) − T q · ∇ T ≥ . (13)The rate in time of the Helmholtz free energy ψ , eq. (10), can be expressed in terms ofits state variables: ˙ ψ = ∂ψ e ∂ (cid:15) e ˙ (cid:15) e + ∂ψ T ∂T ˙ T + ∂ψ pore ∂ (cid:98) ρ ( (cid:98) ρ ) · . (14)The additive decomposition of the strain into an elastic and a (visco-)plastic part,eq.(2), can be inserted into eq.(14) and the result substituted into eq.(13) to yield − (cid:18) ρ ∂ψ T ∂T + ρη (cid:19) ˙ T − T q · ∇ T (cid:124) (cid:123)(cid:122) (cid:125) thermal dissipation + (cid:18) − ρ ∂ψ e ∂ (cid:15) e + σ (cid:19) · ˙ (cid:15) e (cid:124) (cid:123)(cid:122) (cid:125) elastic dissipation + (cid:18) ρ (cid:98) ρ ∂ψ pore ∂ (cid:98) ρ I + σ (cid:19) · ˙ (cid:15) p (cid:124) (cid:123)(cid:122) (cid:125) (visco − )plastic dissipation ≥ . (15)Setting the thermal dissipation to be independent of ˙ T yields equation (9) , whileimposing the vanishing of the elastic dissipation for every ˙ (cid:15) e in equation (15) providesequation (9) , so that using the following simple expression for the potential (borrowedfrom the linear theory) of the elastic part of the free energy density, ρψ e = 12 λ (tr (cid:15) e ) + µ (cid:15) e · (cid:15) e − K b α ( T − T ) tr (cid:15) e , (16)the stress is given by the usual isotropic thermoelastic relation σ = C [ (cid:15) e ] − K b α ( T − T ) I , (17)where K b is the elastic bulk modulus, α is the thermal expansion coefficient, T a referencetemperature, and the fourth-order elastic tensor C is C = λ I ⊗ I + 2 µ S , (18)in which S is the fourth-order symmetrizer and λ and µ are the Lamé elastic moduli.Equation 18 represents a strong assumption, which is justified in the present context,because the elastic strain and rotation are usually small during thermoplastic pressingand sintering of ceramics.As a conclusion, the dissipation inequality eq. (15) reduces to − T q · ∇ T + (cid:18) ρ (cid:98) ρ ∂ψ pore ∂ (cid:98) ρ I + σ (cid:19) · ˙ (cid:15) p ≥ . (19)4 .1 Effective stress for sintering and dissipation An inspection of eq. (19) and consideration of eq. (9) reveals that the thermodynamicdual force to the plastic strain rate is not the stress, but an effective stress defined as [20,21] ˆ σ = σ − σ s I , (20)where σ s = (cid:98) ρR = − (cid:98) ρρ ∂ψ∂ (cid:98) ρ (21)is the so-called ‘sintering stress’ (also known as ‘Laplace pressure’ [31]). Eventually, thedissipation inequality, eq. (19), can be rewritten as − T q · ∇ T + ˆ σ · ˙ (cid:15) p ≥ , (22)which highlights the fact that the inequality is always a-priori satisfied, when the Fourierlaw of heat conduction is assumed q = − k ∇ T, (23)(in which k > is the thermal conductivity), together with the normality rule for ˙ (cid:15) p andconvexity of the yield function, both in the ˆ σ -space. Following [20], the Helmholtz free energy related to the porosity, ψ pore , can be assumedto be a linear function of the surface tension γ s relative to the area of the pore A pore ( (cid:98) ρ ) as follows ψ pore = γ s A pore ( (cid:98) ρ ) ρV = γ s A pore ( (cid:98) ρ ) ρ fd V solid . (24)Assuming spherical pores of radius r in a cubic unit cell, it follows that A pore = 4 πr andassuming incompressibility of the solid phase of volume V solid , r can be expressed as r = (cid:18) V solid − (cid:98) ρ )4 π (cid:98) ρ (cid:19) / , (25)so that the pore potential, eq. (24), becomes ψ pore = γ s ρ fd V / solid π (cid:18) − (cid:98) ρ )4 π (cid:98) ρ (cid:19) / . (26)Finally, a derivative of eq. (26) with respect to (cid:98) ρ yields the sintering stress as σ s = 8 π (cid:18) π (cid:19) (cid:124) (cid:123)(cid:122) (cid:125) ≈ . γ s V / solid (cid:18) (cid:98) ρ − (cid:98) ρ (cid:19) / , (27)where V / solid can be regarded as a length scale, equal to the particle diameter.5 Visco-Plasticity
The existence of a yield function is assumed, depending on the effective stress ˆ σ , whichis the thermodynamic dual of the plastic strain rate, and on the internal variable (cid:98) ρ F = F ( ˆ σ , (cid:98) ρ, T ) . (28)The associative flow rule is assumed [29, 32] involving the non-negative plastic multiplier ˙ λ ˙ (cid:15) p = ˙ λ Q , (29)where the unit yield function gradient is Q = ∂ F (ˆ σ , (cid:98) ρ,T ) ∂ ˆ σ (cid:13)(cid:13)(cid:13) ∂ F (ˆ σ , (cid:98) ρ,T ) ∂ ˆ σ (cid:13)(cid:13)(cid:13) . (30)For rate-dependent problems the coefficient ˙ λ can be replaced with an overstress function,as described in [33]. A possible choice for the overstress function is the yield functiondivided by the viscosity ( η v ). The expression for the strain rate becomes: ˙ (cid:15) p = 1 η v (cid:104) F (cid:105) Q . (31) The Bigoni-Piccolroaz (‘BP’ in the following) yield function [25, 34] is defined by sevenparameters and displays the necessary flexibility to describe a wide range of materialbehaviours and, in particular, can excellently fit the material behavior during the differentstates of sintering. Therefore, the sintering process can be simulated from the beginning,when the material is in a granular form, and then during the firing, up to the finishedceramic piece. The yield function is described by F ( σ , M, m, p c , c, α bp , Θ c , γ bp , β bp ) = F ( p, M, m, p c , c, α bp ) + q g (Θ c , γ bp , β bp ) , (32)where p = − tr σ / , q = (cid:114)
32 dev σ · dev σ , (33)( dev denotes the deviator part) and F ( p, M, m, p c , c, α bp ) = − M p c (cid:113) [Φ bp − (Φ bp ) m ][2(1 − α bp )Φ bp + α bp ] , Φ bp = p + cp c + c ,g (Θ c , β bp , γ bp ) = cos (cid:20) β bp π −
13 cos − ( γ bp cos (3Θ c )) (cid:21) , (34)in which the Lode angle Θ c is defined as Θ c = 13 cos − (cid:18) σ ) q (cid:19) . (35)6arameters α bp and M can be used to adjust the yield function to materials withinternal friction. Parameters β bp and γ bp determine the shape of the yield surface in thedeviatoric plane, so that for both simplicity and lack of experimental data they are setto be zero, which means that Lode angle dependence is excluded (note that for triaxialcompression the Lode angle Θ c is equal to π/ , while it is zero for triaxial extension). Thecohesion c is the hydrostatic yield strength in extension whereas p c is the yield strengthin hydrostatic compression. To model the compaction behavior of a green body during the cold forming process, thefollowing micromechanic assumptions are introduced, which allows the determination ofthe hardening parameters p c , c , and M (dictating the shape of the yield function), asfunctions of the internal variable (cid:98) ρ . It should be noted that the determination of theevolution of parameter p c with the inelastic density (cid:98) ρ is equivalent to the determinationof the theoretical compaction curve for a granulate material. This material is assumed toconsist of elastic perfectly plastic particles, idealized as cylindrical (initially with circularcross section) in a two-dimensional approximation and obeying the Tresca yield criterion,that are equal in size and ordered in a centered cubic geometry, a disposition maintainedfixed during sintering. The results of the two-dimensional approach will be tested againstnumerical simulations for a cubic disposition of perfectly-plastic and initially sphericalparticles, obeying the von Mises criterion. The particles thus yield all simultaneously, sothat to predict the collapse load for this ensemble of particles, the upper bound theoremof limit analysis is applied. A plane strain problem is considered, so that the ceramic particles are assumed in the formof circular cylinders, a geometry which may seem unrealistic, but has been proven to yieldquite reasonable results [12, 35]. Employing the upper bound theorem of limit analysis[36], a collapse mechanism has to be assumed, which, even if it does not correspond toreality, it provides a simple evaluation of the dissipation that occurs within the body.Because the problem is reduced to two-dimensions, the particles are treated as planarfigures, initially circular and later deforming and developing contacts through lines.A collapse kinematics is assumed, in which rigid parts are subject to a pure rigid-bodytranslation, while deformable parts are subject to a pure compressive strain, so that thelatter (denoted by | (cid:15) | in Figure 1) are located on the surface of the particle in contact withthe neighboring particles and the former are localized at the circular corners. The plasticmechanism has two axes of symmetry (vertical and horizontal) so that only a quarter ofthe particle is sketched and the two deformable parts suffer the same strain and dissipatethe same power.The plastic dissipation is generated by the strain rate produced under uniaxial stressin two equally deforming blocks (denoted by | (cid:15) | in Figure 1) and their sliding against: (i.)a square rigid block located at the center of the particle and (ii.) other two rigid blockshaving the shape of a quarter of circle. 7or an equibiaxial (vertical and horizontal) compression load P , which simulates me-chanical compaction, the ‘external’ rate of energy dissipation ˙ W can be written as ˙ W = 2 P v b , (36)where v b is the velocity of the boundary in contact with the neighboring particles. | ε | | ε | v b v b v l v l a h | ε | | ε |a hh=R | ε |a Figure 1: Assuming that the granules behave as rigid-perfectly-plastic materials, theisostatic compaction curve can be calculated using the upper bound technique of limitanalysis. A representative volume element of an idealized 2D granular arrangement isshown in the figure at various stages of the compaction. As the boundaries of the RVEare displaced with the velocities v b , the RVE shrinks. A section view of the assumedcollapse mechanism, for a quarter of the circular particle, is reported in the lower line,where the dashed lines sketch the particle boundaries subject to a displacement rate.The rate of internal power (denoted by D ) is the sum of three dissipation sources, D comp , D slb , and D tr , all related to the strain rate of the blocks denoted with the label | (cid:15) | in Figure 1. These blocks are subject to a linear velocity field preserving incompressibility,so that, assuming a x - x reference system, with x parallel to v b v = v l a x , v = v b h x , (37)where v l is the ‘lateral’ velocity of the block, so that incompressibility requires v b a = − v l h .Therefore, the three sources of dissipation rate can be calculated as follows. • Internal dissipation due to the strain rate of the rectangular blocks | (cid:15) | D comp = 4 kv b a, (38)where k is the limit yield stress under shear (equal to 1/2 of the limit uniaxialstress). 8 Dissipation due to the sliding of the blocks | (cid:15) | against the rigid square block at thecenter of the particle and against the quarter of circle rigid element D slb = 2 k (cid:18)(cid:90) a v dx + (cid:90) h v dx (cid:19) = v l ak + v b hk. (39) • Dissipation due to rigid translation of the quarter of circle rigid element, occurringwith sliding against one of the deforming blocks | (cid:15) | D tr = 2 v l hk. (40)In conclusion, summing up all the contributions from different dissipation sources, theinternal expended power can be written as D = kv b (cid:0) a + a /h + h (cid:1) . (41)Eventually, the limit load P can be expressed as a function of the yielding shear stress k and of the current geometry (specified by a and h ) as P = k (cid:18) a + 12 a h + 12 h (cid:19) , (42)corresponding to the limit pressure that the particle characterized by the dimensions a and h can sustain.The width of the compressed part a and its height h can be expressed as functionsof the relative density and the current side R = a + h of the unit cell (initially equal tothe radius R ), so that, keeping into account incompressibility, the following expressionis found (cid:98) ρ = A A UC = 2 ha + a + πh R = 2 h ( R − h ) + ( R − a ) + πh R , (43)where A UC denotes the area of the unit cell containing the particle of area A . Anexpression for R can also be found using the incompressibility constraint, and thus thefact that the initial area of the two-dimensional particle remains constant: (cid:98) ρ = πR R , ⇒ R = R (cid:114) π (cid:98) ρ . (44)Equation (43) can be solved for h , so that using eq. 44 yields h = (cid:115) πR (1 − (cid:98) ρ ) (cid:98) ρ (4 − π ) (45)and therefore the contact length a results as a = R − h. (46)The hydrostatic yield stress p c is given by the force per unit length P divided by R ,so that the following representation of the compaction curve is obtained:9 c = k √ π (cid:16) − π − (cid:98) ρ ) (cid:113) (cid:98) ρ − π − + 8 (cid:98) ρ (cid:17) (cid:98) ρ ( (cid:98) ρ − . (47)Eq. (47) describes the compaction curve of a granular material and is depicted in Figure3. Formally, Eq. (47) merely represents an upper bound, calculated for a plane strainsituation, which may be believed to be far from reality, so that an assessment of thecapability of eq. (47) to correctly describe a compaction curve is provided through acomparison to a Finite Element simulation involving a three-dimensional distribution ofinitially spherical particles. The simulation was performed with perfectly-plastic particles,ordered in a simple cubic geometry. Initially, the contact between spheres occurs at apoint, but due to the strain, the contact boundary becomes circular and towards thevery end of the compaction tends to become square. This periodic particle arrangementcan be modelled using only one eighth of a sphere, due to symmetry. A FEM modelwas built in Abaqus 6.13, where displacement was imposed through analytical rigid bodyelements, sketched as colored planes in Figure 2, where a unit cell containing a deformedparticle is shown. The material obeys von Mises plasticity, together with a neo-Hookeandescription for the elastic part, with a ratio of Young modulus over yield stress of ,to approximate the rigid-plastic limit. It can be seen that the solution of the limit analysisis, at low densities (about 0.8 and less), remarkably different from the numerical solutionof the 3D problem, which is a direct consequence of the fact that the initial density for the2D problem is different than that pertinent to the 3D arrangement. In fact, the relativedensity in plane strain is about 78%, while this density decreases to about 52%, in a cubicdisposition of spherical particles. Moreover, a typical ceramic powder (as that later usedfor experiments, section 7.3) has an initial density of about 38%. A very simple way forcorrecting the discrepancy between the plane strain upper bound estimate and the valuestypical of ceramic powders is to introduce a correction factor ζ multiplying h , so thatthe initial value of p c becomes correct even for initial densities of 38% and therefore thefollowing ‘modified’ value for h is proposed h mod = ζ R (cid:115) π (1 − (cid:98) ρ ) (cid:98) ρ (4 − π ) , (48)where ζ = 2 . provides the best fit to the material considered in the experiments. Asubstitution of eq. (48) into eq. (36) yields an analytical expression for the compactionbehaviour p c ( (cid:98) ρ, σ m ) = σ m −
530 + (619 + 25 π − (cid:98) ρ ) (cid:113) (cid:98) ρ − π − + 530 (cid:98) ρ √ (cid:98) ρ − . (49)where σ m = k √ is the uniaxial stress for yielding. Note that eq. (49) provides ananalytical description for the compaction curve of a granulate, which will be referred asthe ‘modified limit analysis solution’.Figure 4 shows the Gurson [37] and Helle/Fleck [12, 38] models compared with thepreviously-described FE simulation of a three-dimensional cubic disposition of initiallyspherical particles (‘FEM’ in the figure) and the ‘modified limit analysis solution’ (‘MLA’in the figure). The modified limit analysis solution lies below the curves correspondingto the Gurson and Helle/Fleck models, at high porosity. It needs to be mentioned thatthe Gurson model was originally developed for materials with high density, where the10igure 2: A finite element simulation of a cubic disposition of initially sherical particles,where, due to symmetry, only one eight of one particle is considered. The particle ispressed into a cubic shape by three analytical contact surfaces (shown as translucent redsurfaces) that move towards the center of the cell, with imposed displacements. Incom-pressible hybrid elements are used together with a neo-Hookean material model for theelastic part, while the plastic behaviour is modelled with von-Mises perfect plasticity.11 .0 0.2 0.4 0.6 0.8 1.00246810 Figure 3: Comparison between the two compaction curves obtained from a two-dimensional application of the upper bound theorem of limit analysis and a three-dimensional finite element simulation of the unit cell shown in Figure 2, representative ofthe behavior of a cubic disposition of initially spherical particles. A deformed particle atdifferent levels of density is shown. 12 .0 0.2 0.4 0.6 0.8 1.00123456
Figure 4: Comparison of the Gurson ([37]) and Helle/Fleck [12, 38] models (for thecompaction curve of a porous material) with the modified upper bound solution, ‘MLA’,eq. (49), and the three-dimensional finite element simulation, ‘FEM’, of the unit cellshown in Figure 2. 13odified upper bound solution and the numerical simulation almost coincide, whereasthe model presented by Helle [38] (and used also by Fleck [12]) was intended for cases ofhigh porosities, where predictions of this model are consistent with the modified upperbound and the numerical simulation.
The values for p c (under compressive load) have been determined as a function of (cid:98) ρ .The cohesion strength c can also be evaluated as a function of the internal variable (cid:98) ρ from the following simple micromechanical model. Simulations performed with the unitcell sketched in a deformed state in Fig. 2 show that the contact area between initiallyspherical grains in a cubic geometry is a complex function of the applied pressure, sothat an initially circular contact area evolves towards a square contact shape in the limitof a fully dense material. Therefore, the estimation of the contact areas between grainsin a real powder subject to cold pressing is awkward. Nevertheless, following [38], anapproximation for the normalized (and dimensionless) contact area A c between grainscan be assumed as A c = 4 π (cid:98) ρ − (cid:98) ρ − (cid:98) ρ , (50)so that the stress needed to produce yielding under hydrostatic tension (which is thedefinition of the cohesion c ) is given by the elementary formula c = σ m A c , (51)where σ m is the yielding stress in uniaxial tension of the grains. Note that the normalizedcontact area goes to unity as (cid:98) ρ approaches one. Using again relation (50), the parameter M of the BP yield surface can be expressed as afunction of the internal variable (cid:98) ρ . In particular, it is assumed that at failure under pureshear ( p = 0 so that σ = dev σ ) the deviatoric component of the stress tensor becomes,according to the von Mises criterion, σ m A c / , so that the stress invariant q at failure(denoted by q s ) reduces to σ m A c . Therefore, the following relation is obtained (assuming β = γ = 0 ) for the parameter MM = √ σ m A c p c (cid:113) [ cp c + c − ( cp c + c ) m ][2(1 − α bp )( cp c + c ) + α bp ] . (52) The consolidation process of a ceramic body is activated at high temperatures, wherethe material exhibits a strong thermal softening. Also, the viscosity changes drasticallywith temperature, so that viscoplastic processes occur at much higher velocities, whenthe temperature is high. 14 .1 Temperature evolution
In the framework proposed by Simo and Miehe [22], the temperature evolution is governedby the Fourier law of heat conduction, together with a maximum dissipation postulate.In the absence of internal heat sources, the first law of thermodynamics is ˙ T = 1 ρc h (cid:16) div ( k ∇ T ) + H + D mech (cid:17) , (53)where c h is the specific heat at constant values of strain and internal variables, eq. (12), k the thermal conductivity, eq. (23), D mech is the mechanical dissipation power, D mech = σ · ˙ (cid:15) p + R ( (cid:98) ρ ) · = ˆ σ · ˙ (cid:15) p , (54)and H is the non-dissipative elastic-plastic heating H = T (cid:20) ∂ ψ∂T ∂ (cid:15) e · ˙ (cid:15) e + ∂ ψ∂T ∂ (cid:98) ρ ( (cid:98) ρ ) · (cid:21) . (55)The first term inside the parenthesis describes the piezocaloric effect. The second part isa coupling term between the temperature and the internal variable (relative density). Itis common for plasticity models [22] to replace the mechanical dissipation term with D mech = X ˆ σ · ˙ (cid:15) , (56)where X is sometimes called Quinney-Taylor coefficient and lies between . and . .This factor was introduced as it was found in experiments [39] that the heating associatedwith the plastic flow is roughly 0.9 times the rate of plastic working. If the plastic strainbecomes very large, the Quinney-Taylor factor approaching unity, so that eq. (56) reducesto eq. (54). However, the mechanical dissipation plays only a negligible role in hightemperature applications. In the case of ceramic sintering, a substantial amount of heatis transferred to the system from outside (for instance, the ceramics is left 30 min inthe oven at about 1200 ◦ C), whereas the plastic deformation of the body occurs at arelatively slow pace. This heat transfer is far larger than the heating introduced by thepiezocaloric effect or the plastic heating. It is therefore reasonable to neglect these terms(the mechanical dissipation and the non-dissipative elastoplastic heating) and assume X = 0 , which leads to a system where temperature evolution is uncoupled from plasticflow. Then, eq. (53) reduces to: ˙ T = 1 ρc h div ( k ∇ T ) . (57) At high temperature, a significant thermal softening effect occurs, modelled through avariation of the parameter p c , which is assumed to follow the rule p c ( (cid:98) ρ, σ m , T ) = f T ( T ) p c ( (cid:98) ρ, σ m ) , (58)with f T ( T ) being a steadily decreasing function of temperature, equal to C T at areference temperature of ◦ C . For high temperatures, over about ◦ C and slow defor-mation rates, the material exhibits creep as dominant deformation mechanism [40, 41],and its strength is assumed to become negligible, so that f T ( T ) = (cid:28) − TT C (cid:29) b + C T , (59)15here the Macaulay bracket has been used and T C = 800 ◦ C and C T = 0 . , with b = 0 . . The sintering process is usually accompanied by grain coarsening [42], which has someeffect on the process kinetics and therefore has to be considered in the modelling. Thegrain growth occurring during sintering is assumed to follow an exponential law thatmodels the evolution of the average grain size, as described in [3] and [43], namely ˙ R = γ b M gc R , M gc = M gc exp (cid:18) − Q gc R g T (cid:19) , (60)where R g is the universal gas constant (8.314 J/(mol K)), the constant M gc is taken from[3] to be equal to 2.25 m s/kg and Q gc is the activation energy for grain coarsening. Thevalue for the grain boundary energy γ b is also taken from [3] and is equal to 1.10 J/m .The values of the activation energies, to be used for the model, can be obtained by fittingthem to experimental results, with a procedure described in section 7.1. The shrinkage of the ceramics under thermal load during sintering is essentially a time-dependent process, so that the previously developed time-independent model has to beenhanced to describe rate-dependent effects, thus introducing viscosity. There are sev-eral possibilities to introduce viscous behaviour, for instance using the Coble [44] or theNabarro-Herring [45] creep model. Olevsky [6] introduces a model based on a dissipationpotential that can be derived from a strain energy function. The viscous models proposedin [45, 46] are compared with experimental results in [47, 48] and show relatively good fit.Thus, a viscous model is also used in the present study. The dependence of the viscosity η v on the temperature, density and grain size is assumed to be multiplicative as η v = η v · η v ( R ) · η v ( T ) , (61)where η v is a constant and η v and η v are respectively functions of the radius R of theparticles and of the temperature T . In particular, introducing the initial radius R of theparticles, the function η v is assumed in the form η v ( R ) = (cid:18) RR (cid:19) w , (62)where the exponent w is equal to 3 for grain boundary diffusion as the main type ofinter-particle diffusion mechanism [45]. Finally, the temperature dependence in the law(61) is assumed to follow an Arrhenius type law (compare with [48]): η v ( T ) = e QERgT , (63)where Q E is the activation energy for the viscosity and R g the universal gas constant(8.314 J/(mol K )).In conclusion, the complete function governing the viscosity is η v = η v (cid:18) RR (cid:19) w e QERgT , (64)16hich is used for the viscoplastic flow rule, eq. (31). The values of the constants werefitted against a sintering curve obtained from one dilatometric test, as described in section7.1. The material model developed in the previous Sections was implemented in Abaqus 6.13(Dassault Systèmes) through a UMAT routine created with AceGen [49]. An implicitscheme was used for the yield function, as explained in [50], to overcome difficultieslinked to non-convexity, see [51, 52]. The viscoplastic deformation is evaluated using areturn-mapping algorithm (as described in [22, 53, 54]), which was implemented takingadvantage of the Automatic Differentiation technique [50, 55], to obtain the derivatives(first and second) needed in the algorithm. A consistent tangent matrix for the localNewton method, usually employed in plasticity, was obtained (as described in [23, 54]),without explicitly calculating cumbersome derivatives. The viscosity is implemented usinga Perzyna [33] type approach, as described in Section 3. As the model uses logarithmicstrains, and finite deformation, the ‘Nlgeom’ option of Abaqus was used, which provideslogarithmic strains as the input for the user material routine and uses finite deformationkinematics. The calibration of the yield function, using nominal stress, is described insection 7. For the solution of the global scheme, a separated Newton algorithm was used,in which the coupling terms of the global jacobian matrix (the entries that are a result ofdifferentiation with respect to both displacement and temperature) are set to be zero. Thiscan be justified because the problem is weakly coupled (see Sec. 5.1), and, in addition,the computational cost results strongly reduced. The Newton technique was modifiedusing a line search method, to improve convergence. Note that the general accuracy ofthe solution remains unchanged by using the separated solution scheme, while only therate of convergence is slightly affected.The experimental setup, described in Section 7.3, was modelled in Abaqus, where themold and the stamp were modelled as analytical rigid bodies and the displacement wasimposed on the stamp. To model the furnace environment, the temperature curve wasimposed on the boundaries of the green body, obtained in an initial step through themodelling of the cold powder compaction.
Experiments have been performed to calibrate and subsequently validate the materialmodel. In particular, the calibration of the viscosity parameters of the model has beenconducted with a dilatometric test, while the yield function parameters have been eval-uated on the basis of available experiments [56]. Finally, validation of the model perfor-mance was obtained through a comparison between model predictions and the shrinkingof a specially designed ceramic piece, measured during sintering.
The sintering curve can be obtained using a dilatometer, a standard equipment in labo-ratories dealing with ceramics. Using a special dilatometer and through a long series ofexperiments, a method has been proposed for experimentally determining the viscosities17nd activation energies for pure alumina [47, 48]. Instead of this complex procedure,the viscosity parameters of the proposed model were estimated through a comparisonbetween results from dilatometer tests and numerical simulations in which only a singleelement was used. The temperature was assigned as a constraint on all nodes, followingthe temperature curve of the dilatometer. This procedure can be justified because of thevery small dimensions of the specimen, which does not possess a large thermal capacity.The values for the viscosity constant η v and for the activation energies Q E and Q gc werefound by running simulations with different values and iterating towards an optimal fit.This simplified approach was chosen because experimental values for the activation ener-gies of traditional ceramics are not available. It was assumed that the activation energiesfor grain coarsening and viscosity are equal, an assumption based on the fact that thesevalues have been found very similar for ceramic materials (but different from the powderaddressed in the present study [48]). Experimental data for mechanical properties arerare, as their determination is not straightforward [57] and often not of interest for theinvolved industries. However, a yield strength ranging between MPa and
MPa isprovided in [58], so that σ m = 150 MPa has been assumed. Aluminum silicate spray driedpowder (Sacmi I20087) was used, with a theoretical density (when pores are completelyabsent) of about 2.375 g/cm . To obtain the sintering curve from the powder, a dilatome-ter (TAinstruments DIL-831) was used at a heating rate of 30 ◦ C/min, equal to the heatingrate employed in industrial furnaces. Data in Figure 5 have been plotted after the straindue to thermal expansion was subtracted, to yield a graph that is corrected for thermalexpansion, which was measured in a previous trial, with a completely sintered piece. Thebest fit was found for the viscosity constant η v = 1 · − MPa · s and activation energies Q E = Q gc =340 kJ/mol. The BP yield function, eq. (32), requires the determination of seven parameters, but usingthe previously developed micromechanical modelling, the number of unknown materialconstants can be reduced to four. During powder compaction and sintering, in most ofthe cases the body is under a compressive load. Therefore, the lode angle-dependence isof minor importance and can therefore be neglected, so that parameters β bp and γ bp areset to be zero. The parameters α bp and M determine the shape of the yield surface in the p − q plane. The value m = 4 . has been deduced from [56], while α bp is chosen to beequal to the unit, as in the simple case of the modified Cam-Clay model [25]. The size ofthe yield surface in the p − q plane is determined by p c , c and M , which are functions ofthe relative density. The yield surface evolves as the density grows, so that this evolutionis shown in Figure 6. 18 ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■■■■■■■■■■■■■■■■■■ - - - - - Figure 5: Activation energies and the viscosity constants have been determined withrepeated numerical simulations of the experimental results obtained from a dilatometertest, until a good fit has been found (illustrated in the figure).19
100 0 100 200050100150200
Figure 6: The BP yield surface evolution with respect to the relative density, for a porcel-lain stoneware ceramic (aluminum silicate spray dried powder). The parameters definingthe yield surface are α bp = 1 , m = 4 . , γ bp = 0 , β bp = 0 , σ m = .2 Material parameters used for the examples In summary, the values of the parameters used in the simulations, defining the materialbehaviour and the contact conditions at the powder/mold interface, where a Coulombfriction law was assumed, are listed in Table 7.2.Material Parameters Values R : Initial particle radius 11.24 · − m γ s : Surface energy 1.10 J/m Q E : Viscous activation energy 354 kJ/molM gc : Grain boundary mobility coefficient 2.25 m s/kgQ gc : Grain boundary mobility activationEnergy 354 kJ/mol E : Young’s modulus 5000 MPa X : Taylor Quinney Coefficient 0 σ m : Compressive strengthof the fully dense material 150MPa w : Exponent in the viscosity law 3 R g : Gas constant 8.314 J/(mol K) T C : Temperature constantfor thermal softening law 800 ◦ C T C : Constant for thermal softening law 0.0001 b : Constant for thermal softening law 0.9 ρ : inital density 0.38 η v : viscosity constant − MPa · s m : BP Parameter 4.38 β : BP Parameter 0 γ : BP Parameter 0 α : BP Parameter 1 µ : Die/stamp Coulomb friction coefficient 0.421 .3 Experimental validation Validation of the model has been obtained by referring to the forming and subsequentsintering of a special green piece, namely, a profiled tile. It is possible to introduce somecomplexity to a floor tile by using a special tool with varying heights (as sketched inFigure 7), thus obtaining a green with zones of different height and therefore density.Before sintering, this profiled geometry has to be pressed from powder, using a tool thatwas manufactured for this specific purpose (Fig. 7, lower part). The dimensions beforepressing were 330mm × ◦ C, 1150 ◦ C and 1200 ◦ C, byadopting the heating cycle shown in Figure 8, using an oven available at the Sacmi Labs.The density variation before and after sintering was measured using the purpose-made X-ray Line Scanner (of the ‘Sacmi Continua+’ line, manufactured by Microtec Srl,Bressanone, Italy, Fig. 9) at the Sacmi laboratory. While the thickness changes abruptly,creating a tile with three rather distinct heights (9.8mm, 10.2mm and 10.7mm), thedensity varies less abruptly, as can be seen in Figure 11a. The density distribution wasmeasured on pairs of nominally identical tiles, both subject to the same treatment. Thedensities in the pairs (before and after firing) was found so similar (the experimental datawere almost superimposed), so that only one experiment for each pair is reported below.22igure 7: Upper part: Photo of the green body, with zones of different height and thereforedensity, used for sintering and density measurements. The green was formed with a tool,which was manufactured with three different heights, sketched in the lower part of thefigure (not true to scale). ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆
Figure 8: Sintering was obtained by moving the green through a continuous oven (a‘Sacmi Forni S.pA. EUP 130 was used) across different temperature zones, so that thetime-temperature curve shown above is applied.23igure 9: Density measurements were performed with an X-Ray scanner from the ‘SacmiCONTINUA+’ line, shown in the photo. Note that a profiled tile is entering the scanner.24 .4 Simulation of the forming and sintering of the ceramic platewith different densities
The stamp and mould for powder forming were modelled as rigid bodies in contact withCoulomb friction (coefficient equal to 0.4). The powder was initally considered of uniformheight (22mm) and homogeneous relative density (0.38), as in the experimental setup.The stamp is then displaced by ∆ H = 12.6 mm, pressing the body in the desired shape.Powder pressing is usally modelled as a rate-independent process [15], so that for thepowder pressing part of the simulation, the viscosity was set to a constant and low value,to come close to the limit of rate-independent plasticity. For the simulation of sintering,the viscosity model reported in section 5.4 was used and the temperature curve of theoven was prescribed on the boundaries of the ceramic body. The simulated geometry isshown in Figure 10 (not to scale), where the undeformed piece (a rectangle, unmeshed)and the deformed mesh after pressing and subsequent mold release are reported. TheFigure 10: The mesh of the modelled powder in its initial state (a rectangle, not meshed)and after pressing and subsequent release from the mold (with mesh). Just one half thepiece was simulated, due to symmetry. The figure is scaled by 200% in height direction,to make the contour change more visible.figure clearly shows the large strain suffered by the powder during compaction and thepresence of a modest spring back effect, which is limited in the lateral direction, becauseof the high aspect ratio of the sample. In particular, the tile increases in height afterrelease, but decreases a little bit in width.Measurements of density and thickness are shown and compared to the simulation inFigs. 11 and 12. These figures prove the validity of the model, which is capable ofreproducing the entire compaction and sintering processes of a ceramic granulate with anexcellent precision, so that the percent errors (reported in the lower parts of the graphs)are below 4% at maximum and for most of the part below 2% .Finally, a photograph of the sintered ceramic piece is shown in Figure 13 on the right,while the deformed mesh obtained after the simulated compaction and sintering processisshown on the left. The contour of the photo shown in the figure has been marked in blue,so that it can be observed that the distorsion of the boundary induced by sintering is wellcaptured by the model, another demonstration of the excellent predictive capabilities ofthe mechanical model. 25 ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ (a) Density variation ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ (b) Thickness variation Figure 11: Density and thickness variation of the green after pressing, measured by an X-ray scanner, and compared to the model prediction through numerical simulation. Percenterrors between results from the simulation and the experiment are reported in the lowerparts of the graphs. 26 ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ (a) Density variation ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ (b) Thickness variation
Figure 12: Density and thickness variation of the ceramic piece after firing at 1100 ◦ C(upper part), 1150 ◦ C (central part) and 1200 ◦ C (lower part) measured by an X-ray scan,and compared to the model prediction through numerical simulation. Percent errorsbetween results from the simulation and the experiment are reported in the lower partsof the graphs. 27igure 13: The simulation of the formed and fired ceramic piece (on the left) compared to aphoto of the real ceramic piece, sintered at 1200 ◦ C (on the right, with the contour markedin blue). The qualitative trend of the distorsion of the boundary is well reproduced bythe simulation, which also correctly captures also the pronounced shrinkage at the middleof the ceramic piece. 28
Discussion and conclusions
It has been shown that a thermomechanical model can be formulated, implemented andcalibrated to provide a computational tool for the simulation of the entire process ofceramic production, starting from the cold pressing of a granulate and ending with thesubsequent non-isothermal firing and sintering.Compared with the results of ad hoc performed experiments, the model predictionsprovide an accurate description of the density distribution and the shape distorsion suf-fered by the piece during the production process.Although the developed model is based on several simplificative assumptions and thelack of experimental data (inherent to the extreme conditions to which a ceramic pieceis exposed) has precluded a fine calibration of model parameters, it is believed that thepresented results show that mechanical modelling is a valuable alternative to the empiricalprocesses still often in use for the design of pieces in the ceramic industry.
Acknowledgements.
The authors would like to acknowledge invaluable assistance andguidance during experiments provided by Ing. Claudio Ricci (Sacmi, Imola) and Dr. Si-mone Sprio (CNR-Istec, Faenza). The authors are grateful to Professor Rebecca Brannon(University of Utah) for her kind advice, particularly on the treatment of viscoplasticconstitutive equations.Financial support from FP7-PEOPLE-2013-ITN Marie Curie ITN transfer of knowl-edge programme PITN-GA-2013-606878-CERMAT2 (D.K.), from PRIN 2015 ’Multi-scalemechanical models for the design and optimization of micro-structured smart materialsand metamaterials’ 2015LYYXA8-006 (D.B.) and from the ERC advanced grant ERC-2013-ADG-340561-INSTABILITIES (A.P.) is gratefully acknowledged.29 eferences [1] O. Guillon, J. Rödel, and R. K. Bordia. “Effect of green-state processing on thesintering stress and viscosity of alumina compacts”. In:
Journal of the AmericanCeramic Society doi : .[2] A. T. Procopio and A. Zavaliangos. “Simulation of multi-axial compaction of gran-ular media from loose to high relative densities”. In: Journal of the Mechanics andPhysics of Solids doi : .[3] A. Wonisch, T. Kraft, M. Moseler, and H. Riedel. “Effect of different particle sizedistributions on solid-state sintering: A microscopic simulation approach”. In: Jour-nal of the American Ceramic Society doi : .[4] M. Abouaf, J. L. Chenot, G. Raisson, and P. Bauduin. “Finite element simulationof hot isostatic pressing of metal powders”. In: International Journal for NumericalMethods in Engineering doi : .[5] J. M. Duva and P. D. Crow. “The densification of powders by power-law creep duringhot isostatic pressing”. In: Acta Metallurgica et Materialia doi : .[6] E. A. Olevsky. “Theory of sintering: from discrete to continuum”. In: MaterialsScience and Engineering: R: Reports doi : .[7] T. Kraft and H. Riedel. “Numerical simulation of solid state sintering; model andapplication”. In: Journal of the European Ceramic Society doi : .[8] S Stupkiewicz, A. Piccolroaz, and D. Bigoni. “Finite-strain formulation and FE im-plementation of a constitutive model for powder compaction”. In: Computer Methodsin Applied Mechanics and Engineering
283 (2015), pp. 856–880. doi : .[9] M Scot Swan, A. Piccolroaz, and D. Bigoni. “Application of tomographic reconstruc-tion techniques for density analysis of green bodies”. In: Ceramics International doi : .[10] M Penasa, L Argani, D Misseroni, F Dal Corso, M Cova, and A Piccolroaz. “Com-putational modelling and experimental validation of industrial forming processes bycold pressing of aluminum silicate powder”. In: Journal of the European CeramicSociety doi : .[11] N. A. Fleck. “On the cold compaction of powders”. In: Journal of the Mechanics andPhysics of Solids doi : .[12] N. Fleck, L. Kuhn, and R. McMeeking. “Yielding of metal powder bonded by isolatedcontacts”. In: Journal of the Mechanics and Physics of Solids doi : .[13] P. Redanz and N. A. Fleck. “The compaction of a random distribution of metal cylin-ders by the discrete element method”. In: Acta Materialia doi : .3014] R. W. Lewis, D. T. Gethin, X. S. Yang, and R. C. Rowe. “A combined finite-discreteelement method for simulating pharmaceutical powder tableting”. In: InternationalJournal for Numerical Methods in Engineering doi : .[15] A. Piccolroaz, D. Bigoni, and A. Gajo. “An elastoplastic framework for granularmaterials becoming cohesive through mechanical densification. Part I – small strainformulation”. In: European Journal of Mechanics - A/Solids doi : .[16] A. Piccolroaz, D. Bigoni, and A. Gajo. “An elastoplastic framework for granularmaterials becoming cohesive through mechanical densification. Part II – the formu-lation of elastoplastic coupling at large strain”. In: European Journal of Mechanics- A/Solids doi : .[17] A. C. Cocks and I. C. Sinka. “Constitutive modelling of powder compaction - I.Theoretical concepts”. In: Mechanics of Materials doi : .[18] A Salvadori, S Lee, A Gillman, K Matouš, C Shuck, A Mukasyan, M. Beason, I.Gunduz, and S. Son. “Numerical and experimental analysis of the Young’s modulusof cold compacted powder materials”. In: Mechanics of Materials
112 (2017), pp. 56–70. doi : .[19] A Krairi, K Matouš, and A Salvadori. “A poro-viscoplastic constitutive model forcold compacted powders at finite strains”. In: International Journal of Solids andStructures
135 (2018), pp. 289–300. doi : .[20] L. Mähler and K. Runesson. “Constitutive Modeling of Cold Compaction and Sin-tering of Hardmetal”. In: Journal of Engineering Materials and Technology doi : .[21] J Frischkorn and S Reese. “Simulation of powder metallurgical coating by radialaxial ring rolling”. In: Pamm doi :
10 . 1002 / pamm .200810525 .[22] J. C. Simo and C Miehe. “Associative coupled thermoplasticity at finite strains:Formulation, numerical analysis and implementation”. In:
Computer Methods inApplied Mechanics and Engineering doi : .[23] C. Miehe, N. Apel, and M. Lambrecht. “Anisotropic additive plasticity in the loga-rithmic strain space: Modular kinematic formulation and implementation based onincremental minimization principles for standard materials”. In: Computer Meth-ods in Applied Mechanics and Engineering doi : .[24] C. Sansour and W. Wagner. “Viscoplasticity based on additive decomposition of log-arithmic strain and unified constitutive equations - Theoretical and computationalconsiderations with reference to shell applications”. In: Computers and Structures doi : .[25] D. Bigoni and A. Piccolroaz. “Yield criteria for quasibrittle and frictional materials”.In: International Journal of Solids and Structures doi : . arXiv: .3126] S Stupkiewicz, A Piccolroaz, and D Bigoni. “Elastoplastic coupling to model coldceramic powder compaction”. In: Journal of the European Ceramic Society doi : .[27] A. Gajo and D. Bigoni. “Elastoplastic coupling for thermo-elasto-plasticity at hightemperature”. In: Geomechanics for Energy and the Environment doi : .[28] L. P. Argani, D Misseroni, A Piccolroaz, Z Vinco, D Capuani, and D Bigoni.“Plastically-driven variation of elastic stiffness in green bodies during powder com-paction: Part I. Experiments and elastoplastic coupling”. In: Journal of the EuropeanCeramic Society doi : .[29] I. Doghri. Mechanics of Deformable Solids: Linear, Nonlinear, Analytical, and Com-putational Aspects . Berlin: Springer, 2000.[30] G. A. Maugin.
The thermomechanics of plasticity and fracture . Cambridge: Cam-bridge University Press, 1992.[31] L. Galuppi and L. Deseri. “Combined effects of interstitial and Laplace pressure inhot isostatic pressing of cylindrical specimens”. In:
Journal of Mechanics of Materialsand Structures doi : .[32] J. Lubliner. Plasticity Theory . Mineola: Dover, 2008, p. 528.[33] P. Perzyna. “The constitutive equations for rate sensitive plastic materials”. In:
Quarterly of Applied Mathematics
XX.4 (1963), pp. 321–332.[34] A. Piccolroaz and D. Bigoni. “Yield criteria for quasibrittle and frictional materials:A generalization to surfaces with corners”. In:
International Journal of Solids andStructures doi : .[35] L. P. Argani, D. Misseroni, A. Piccolroaz, D. Capuani, and D. Bigoni. “Plastically-driven variation of elastic stiffness in green bodies during powder compaction: PartII. Micromechanical modelling”. In: Journal of the European Ceramic Society doi :
10 . 1016 / j . jeurceramsoc . 2016 . 02 . 012 . arXiv: .[36] W. F. Chen.
Limit analysis and soil plasticity . Ross, 2008.[37] A. L. Gurson. “Continuum Theory of Ductile Rupture by Void Nucleation andGrowth: Part I-Yield Criteria and Flow Rules for Porous Ductile Media”. In:
Jour-nal of Engineering Materials and Technology doi : .[38] A. S. Helle, K. E. Easterling, and M. F. Ashby. “Hot-isostatic pressing diagrams:New developments”. In: Acta Metallurgica doi : .[39] G. I. Taylor and H. Quinney. “The Latent Energy Remaining in a Metal afterCold Working”. In: Proceedings of the Royal Society A: Mathematical, Physical andEngineering Sciences doi : .arXiv: .[40] M. F. Ashby and H. Frost. Deformation-mechanism maps . Cambridge UniversityEngineering Dept., 1981. 3241] V. H. Hammond and D. M. Elzey. “Elevated temperature mechanical properties ofpartially sintered alumina”. In:
Composites Science and Technology doi : .[42] M. N. Rahaman. Sintering of ceramics . CRC Press, 2007.[43] M Hillert. “On the theory of normal and abnormal grain growth”. In:
Acta Metallur-gica doi : . arXiv: arXiv:1011.1669v3 .[44] R. L. Coble. “A Model for Boundary Diffusion Controlled Creep in PolycrystallineMaterials”. In: Journal of Applied Physics doi :
10 .1063/1.1702656 . arXiv: arXiv:1011.1669v3 .[45] M. Rahaman, L. Jonghe, and R. Brook. “Effect of Shear Stress on Sintering”. In:
Journal of the American Ceramic Society doi : .[46] R. K. Bordia and G. W. Scherer. “On constrained sintering-I. Constitutive modelfor a sintering body”. In: Acta Metallurgica doi : .[47] R. Zuo, E. Aulbach, and J. Rödel. “Experimental determination of sintering stressesand sintering viscosities”. In: Acta Materialia doi : .[48] R. Zuo and J. Rödel. “Temperature dependence of constitutive behaviour for solid-state sintering of alumina”. In: Acta Materialia doi : .[49] J. Korelc. “Automation of primal and sensitivity analysis of transient coupled prob-lems”. In: Computational Mechanics doi :
10 . 1007 /s00466-009-0395-2 .[50] S Stupkiewicz, R Denzer, A. Piccolroaz, and D. Bigoni. “Implicit yield functionformulation for granular and rock-like materials”. In:
Computational Mechanics doi : .[51] R. M. Brannon and S Leelavanichkul. “A multi-stage return algorithm for solvingthe classical damage component of constitutive models for rocks, ceramics, and otherrock-like media”. In: International Journal of Fracture doi : .[52] M Penasa, A Piccolroaz, L Argani, and D Bigoni. “Integration algorithms of elasto-plasticity for ceramic powder compaction”. In: Journal of the European CeramicSociety doi : .[53] J. C. Simo and T. J. R. Hughes. Computational Inelasticity . New York: Springer,1998.[54] E. A. de Souza Neto, D Perić, and D. R. J. Owen.
Computational Methods ForPlasticity . West Sussex: John Wiley & Sons, Ltd, 2008.[55] J. Korelc and P. Wriggers.
Automation of Finite Element Methods . Cham: SpringerInternational Publishing, 2016. doi : .3356] F. Bosi, A. Piccolroaz, M. Gei, F. D. Corso, A. Cocquio, and D. Bigoni. “Experi-mental investigation of the elastoplastic response of aluminum silicate spray driedpowder during cold compaction”. In: Journal of the European Ceramic Society doi :
10 . 1016 / j . jeurceramsoc . 2013 . 11 . 037 . arXiv: arXiv:1405.0932v1 .[57] M. Gei, D. Bigoni, and S. Guicciardi. “Failure of silicon nitride under uniaxial com-pression at high temperature”. In: