Numerical Simulation of the Temperature and Stress Field Evolution Applied to the Field Assisted Sintering Technique

The field assisted sintering technique (FAST) is a high-amperage, low-voltage, powder consolidation technique that employs pulsed direct current and uniaxial pressure. Over the past several years, FAST has been successfully used to produce a variety of different materials including metals, composites, and ceramics. In this paper we present a transient finite element model of aluminum oxide sintering that incorporates a coupled electrical, thermal, and mechanical analysis that closely resembles the procedures used in physical experiments. Within this context, we outline the governing equations that pertain to a balanced energy equation and include the effects of thermal and electrical contact forces, radiation, and Joule heating. We couple this with the relevant equations pertaining to mechanical displacements and prescribe the necessary initial and boundary conditions for a complete solution. As part of our transient analysis, we also present our implementation of a proportional integral derivative controller, which (similar to actual experimental conditions) affords the use of a predetermined heating rate conditioned upon a variable voltage. Finally, we discuss implications relating to the temperature and stress fields and suggest possible avenues for improvement.


Introduction
The field assisted sintering technique (FAST) (also known as spark plasma sintering (SPS) or pulsed electric current sintering (PECS)) is a high-amperage, low-voltage, powder consolidation technique that employs pulsed direct current (DC) and uniaxial pressure [1,2].Over the past several years, FAST has been successfully used to produce a variety of different materials including metals [3], composites [4], and ceramics [5,6].In contrast to the more traditional (pressure assisted) methods such as hot pressing (HP) or hot isostatic pressing (HIP), some of the most significant benefits of using FAST include the reduction in processing times, the elimination of sintering aids, and high heating temperatures [7].This helps minimize grain growth and can result in mechanical [8], physical [9], and optical [10] improvements in the material.
Over the past decade, there has been a significant increase in FAST research, as evidenced by the increase in the number of papers on the subject [7].While a significant proportion of this research activity has been strictly experimental, more recently, there has been an increase in the number of experiments performed in conjunction with numerical modeling efforts.This integration has helped overcome some of the inherent challenges associated with performing experimental observations (i.e., the direct measurement of the temperature and stress field within the sample) and led to a better understanding of the cause and effect relationships between the evolving processing parameters (i.e., stress and temperature fields) and the evolution of the material microstructure.In particular, there have been several recent works that have used the finite element method to investigate the SPS temperature distributions [11][12][13][14][15][16][17][18][19].While each of these studies has provided significant contributions to the field, there still remains much room for improvement.Relatively few of these studies, for example, take into account the presence of contact resistances [11,12,14,15], and the majority fail to account for the stress distribution due to uniaxial  loading.This latter omission is particularly relevant, since several recent experimental studies have demonstrated that the inclusion of processing stresses is of vital importance to densification mechanisms [7,20].An additional omission includes a coupled approach to the applied current (or voltage) and the desired heating rate.Exceptions include the works of Antou et al. [16], Munoz and Anselmi-Tamburini [18], and Wang et al. [19], who each presented a fully integrated thermal, electrical, and mechanical approach to the problem.The aforementioned numerical experiments, however, are all limited solely to SPS research and have neglected (to the best of the knowledge of the present authors) any reference to FAST geometries.
In this work we present a FAST, transient, and finite element model that incorporates a coupled electrical, thermal, and mechanical analysis that closely resembles the actual experimental operating conditions.Within this context, we outline the governing equations that pertain to a balanced energy equation and include the effects of contact forces, radiation, and Joule heating.We couple this with the relevant equations pertaining to mechanical displacements and prescribe the necessary initial and boundary conditions for a complete solution.As part of our transient analysis, we also present our implementation of a proportional integral derivative (PID) controller, which (similar to actual experimental conditions) affords the use of a predetermined heating rate conditioned upon a variable voltage.We focus our attention on a single candidate material: fully dense aluminum oxide (Al 2 O 3 ), and we provide experimental validations where relevant.Finally, we discuss implications relating to the temperature and stress fields and suggest possible avenues for model improvement.

Method
2.1.Geometry.Schematics showing both a solid and wireframe representation of the FAST apparatus (FCT HP 25/1, FCT, Germany) used in this work are provided in Figures 1(a) and 1(b), respectively.The cylindrically symmetric apparatus is composed of the sample, a graphite die, two opposing punches, and two steppers consisting of opposing circular frustums.The purpose of the enclosed die is to mitigate lateral shrinkage of the sample during the uniaxial mechanical loading and promote pressure-induced densification only along the axial direction.As indicated in Figure 1(a), each of these tool geometries are composed of graphite (type R7710, SGL Carbon GmbH, Germany) [21].The primary dimensions for each component of the apparatus are shown in Table 1 and Figure 1  that has a diameter of 10 mm and is 31 mm deep.The steppers are conical and are placed on either end of the punches.The steppers contain maximum and minimum diameters of 80 mm and 33.8 mm, respectively, and contain a 20 mm recessed hole for the punch to sit in (with a height of approximately 2.5 mm).Additionally, the steppers contain a 10 mm diameter hole that lines up with the holes in the punch.
The central hole axially penetrating both the punches and the steppers is used to accommodate an optical pyrometer that sits approximately 5 mm from the top edge of the sample.For comparison purposes, this location is referenced in Figure 2 as point 5. Points 1-4 correspond to additional reference points that will be addressed in Section 3.

Material Properties.
The temperature-dependent material properties for both the graphite tooling assembly (type R7710, SGL Carbon GmbH, Germany) [21] and the Al 2 O 3 sample compound used in this study are shown in Table 2.As indicated, these values correspond to temperatures between approximately 273 K and 2500 K. Apart from the temperature, the density dependence may also require special consideration.However, since the density of the graphite over this temperature range varies by less than 5% [21], it was assumed constant and modeled as an isotropic material.In contrast, the Al 2 O 3 is known to have a highly variable density that is dependent on the sintering stage [18].This density dependence corresponds to the various diffusion-related processes which contribute to such densification mechanisms as grain growth, pore shrinkage, and pore annihilation.While various mesoscale modeling efforts, including Tikare et al. [22], have made considerable progress in predicting this density evolution, still, the majority are highly dependent on physical experimentation for their initial and boundary conditions.These limitations are further compounded at larger scales, where more of the apparatus is modeled.For these reasons, this work, like the majority of previous studies [11][12][13][14][15][16][17][18][19], assumes the sample to have a constant density, indicative of the final stages of sintering.

Governing Equations.
The set of equations governing the three-way coupled interaction between the electrical, thermal, and mechanical system response are shown in ( 1)- (8).The four dependent variables (degrees of freedom), consisting of temperature, the radial and axial components of displacement, and the electric potential, are solved through an energy balance approach acting over the entire apparatus.As indicated in (1), the relative heat gain by conduction is governed by Fourier's law, while the remaining heat transfer mechanisms, q conv , q j , q rad , and q cont , corresponding to convection, Joule heating, radiation, and contact resistance, respectively, are assigned to source terms.In addition to the contribution from uniaxial loading, the mechanical stress, shown in ( 6)- (8), is also coupled to the energy equation, through thermal expansion.Additionally, since the majority of the material properties are thermally dependent (see Table 2), this requires an iterative approach to convergence where ρ is the sample density, C p is the sample specific heat at constant pressure, T is the temperature, k is the thermal conductivity, and q j , q rad , q conv , and q cont are source terms corresponding to Joule heating, radiative heat transfer, convective heat transfer, and contributions from thermal and electrical contact resistances, respectively.Since the experiments take place within a vacuum, the heat loss due to convection can be neglected (q conv = 0).Surface to ambient radiation occurs among the exposed surfaces of the die, spacers, and punches and is given by where σ s is the Stefan-Boltzmann constant (σ s = 5.67•10 −8 W/m 2 K 4 ), A s is the area of the exposed surface, ε is the emissivity (ε graphite = 0.8 [11,12]), and T e and T a are the emitting and absorbing surface temperatures, respectively.The term q j represents the rate of Joule heating per unit volume and is given by where J is the current density and E is the electric field.These in turn are computed via solutions of (4), where σ is the electrical conductivity and ϕ is the electric potential.
As shown, the electrostatic condition is imposed to reduce modeling complexity and is in general a valid assumption since the electric potential reaches a state of equilibrium far in advance of the temperature.
It has been demonstrated that thermal and electrical resistances induced at both the horizontal and vertical contacts between the various components of the tool assembly and sample can have a significant influence on the evolution of the temperature and electrical potential fields [11,12].While there exist several different approaches to account for these resistances, this work treats them in accordance with Vanmeensel et al. [12].
where h g and σ g are the thermal and electrical contact conductance coefficients, respectively.Furthermore, it was assumed that the vertical contact coefficients (h g|V and σ g|V ) are proportional to the horizontal coefficients (h g|H and σ g|H ) such that h g|H = 15E3 W/m 2 K, σ g|H = 5E7 Ω −1 m −2 , h g|V = h g|H /6, and σ g|H = σ g|H /6.While these values are only best fit approximations, they correspond to the values selected by several other authors [11,12].Subsequent to the solution to the temperature field, the mechanical displacements are computed from the equations of linear elasticity, under the simplifying assumption that the sample yield stress is not exceeded.These are shown in ( 6)- (8) as where σ is the stress tensor, F is the body force per unit volume, ε is the strain tensor, u is the displacement vector, ν is Poisson's ratio, μ is the Lame coefficient (μ = μ(E, ν)), and E is Young's modulus.
In order to solve the governing equations, the finite element method was utilized.The variational form of (1), after integration by parts of the conduction term and the application of the divergence theorem, yields where ν i is the basis function, Ω and Γ are the elemental volume and its enclosing surface, respectively, and n is the unit normal vector to the surface.The linearized partial differential equation for T in the variational formulation is stated as where M i j is the mass matrix, A i j is the stiffness matrix, and F i is the force vector.Comparing (10) with ( 9) leads to the following: where q = k∇T • n is the heat flux perpendicular to the surface normal.

Mesh and Solution Discretization.
As shown in Figure 2, the geometric symmetry of the SPS apparatus allows for a two-dimensional (2D) axisymmetric modeling approach in cylindrical coordinates.As such, the linearized partial differential equations ( 10)- (11) were simplified to reflect this change of coordinates.The 2D domain was meshed via Delaunay triangulation and was composed of 12,428 threenode triangular elements with an average element quality of 0.95.For each node, four degrees of freedom were assigned and consisted of displacement in the r and z directions, temperature, and electric potential.The total number of degrees of freedom was 101,747.The numerical integration over each triangular element was approximated by the sum over all Gauss-point contributions and solved (for each degree of freedom) with an absolute error tolerance of 0.001 using the multiphysics simulation software COMSOL 4.1.Time stepping was performed using a second-order backward differencing scheme.

Boundary Conditions.
After initializing the temperature over the entire domain, T 0 = 673.15K, the boundary conditions were applied in accordance with Figure 2. As stated previously, since the experiments occurred within a vacuum, the heat losses by convection were neglected.A constant temperature of 673.15K was applied to the fixed (zero displacement) bottom stepper surface (due to refrigerant cooling) along with an electric potential, V = 0 (corresponding to the grounding condition).The sides were assigned ambient to surface radiation boundaries implemented in accordance with (2).The top stepper was also assigned a constant temperature of 673.15 K, along with a time dependent applied load (P) which was assigned in accordance with Figure 2(b).The top stepper was also given a voltage (dependent on a prescribed heating rate) and functioned in accordance with a proportional integral derivative (PID) controller described in greater detail below.

PID Controller. FAST experiments in general follow
a preset heating rate that is input by the operator and is regulated by a time changing current or voltage profile.This allows for a high degree of temperature control and thus better overall control over the evolving sample microstructure.
Ideally, a numerical model should follow a similar control process.In fact, however, the majority of the numerical work done previously was conducted under the assumption of a constant applied current [13,17,23,24].In the present work a PID controller was utilized to achieve a close approximation to the experimental process.In general, a PID controller approximates an error correction between a measured process variable and a desired set point by the calculation and subsequent output of a corrective action, which can adjust the process in an efficient manner and maintain minimal error.[18].The PID controller thus works in an iterative and closed-loop configuration.In this work the applied voltage V (t) is adjusted in accordance with a prescribed heating rate (R).A tracking error e(t), representing the difference between the desired heating rate and the actual heating rate (R act ), is computed for a specified control point (see point 5 of Figure 2).This error signal e(t) is then sent to the PID controller, which computes the new voltage via the derivative and the integral of this error signal according to where K P , K I , and K D represent the proportional, integral, and derivative gains, respectively.The constant values for these gains were chosen for an optimal control response and were assigned values in accordance with the Zeigler-Nichols method [15].

Results and Discussion
The physical experiments were performed by controlling the temperature (as measured at point 5 of Figure 2) according to a predetermined temperature-time profile (see Figure 3(a)).
As indicated, the temperature-time curve assumes an initially constant temperature of 673.15K for the first 192 seconds.Thereafter, a linear increase of 3.24 K/s is maintained until approximately 376 s, at which the temperature is again fixed at 1275 K.A second linear increase is performed from 560 to 857 seconds at a rate of 1.68 K/s and peaks at a temperature of 1773 K where there is an approximate dwell period of approximately 364 s.Thereafter, at 1221 seconds, cooling begins at a linearly decreasing rate of 1.68 K/s and continues until the electrical power is shut off at approximately 1516 seconds.The total duration of the simulated experiments is just over 2200 s (36 min).Unlike many other FAST experiments, the applied uniaxial loading was not constant but varied with time according to Figure 3(b).As indicated, the uniaxial applied force is constant at 5 KN for the first 375 seconds and then linearly increases at a rate of 0.15 KN/s to a maximum of 31 KN.Thereafter, at a time of 1515 seconds the force decreases at a rate of 0.22 KN/s to again match its initial magnitude of 5 KN. Figure 3(a) shows the simulated temperature profile (dashed line) and  the actual experimental temperature (solid line) measured at point 5 of Figure 2. As shown, apart from some minor variation at the peak temperature and at the onset of the electrical power shutoff, the agreement is quite good.While these variations could potentially be mitigated by selecting different combinations of PID gain parameters, due to the excellent agreement in the remaining temperature-time profile this was deemed unnecessary.a drop is observed at time t = 376 s in order to stabilize the controlled temperature at 1300 K.The second heating ramp culminating in a peak voltage (at point 5) of nearly 1.75 Volts is (on average) maintained for the remaining dwell period.The cooling ramp declines almost linearly until the voltage is shut off at 1516 seconds.In contrast, the potential within the aluminum oxide sample (point 1) shows overall a slightly increased voltage culminating in an electric potential of approximately 1.875 V at its peak.
For comparison purposes, the temperature-time profiles at points 2-5 are shown in Figure 4.As indicated, all of the points show similar profiles and vary only with respect to amplitude.This variation occurs primarily between the end of the first heating ramp and the point of electrical shutoff (376 s < t < 1516 s).As shown, the highest temperatures, corresponding to nearly 1800 K, were obtained along the dwell period (857 s < t < 1221 s) and apply to points 5 and 3 located in close proximity along the axial centerline.The minimum temperature within the dwell period, corresponding to approximately 1568 K, was obtained at point 2 and corresponds to the radial interface between the sample and the die.This decrease in temperature may further be illustrated by plotting the temperature along the line connecting points 1 and 2, at a single point of time t = 1000 s, as shown in Figure 4(b).As indicated, the radial temperature difference between the sample centerline and the sample-die interface amounts to approximately 7.5%.As in the work of Munoz and Anselmi-Tamburini [18], this rather significant level of nonhomogeneity may be attributed to the relatively low thermal conductivity of aluminum oxide, as compared to other more conductive sample specimens, which tend to reveal (for equivalent operating conditions) a significantly more uniform temperature distribution.Figures 5(a)-5(d), show temperature contour plots within the sample domain at four different time sequences, corresponding to t = 300 s, t = 600 s, t = 1000 s, and t = 1300 s, respectively.As shown, each of the plots clearly identifies a cold region along the right-hand side corresponding to the sample-die interface.Further, all but Figure 5(d), show the hottest zones at the top and bottom of the sample surface, corresponding to the sample punch interface.Figure 5(d) (t = 1300 s), which corresponds to the cooling ramp, shows, in contrast, a hot zone centrally located along the central axis of the sample.These plots also indicate, with the aid of Figure 5(a) and the plotted isotherms, some useful trends concerning the temperature gradient.As indicated, the highest temperature difference between points 1 and 2 corresponds to the dwell interval (857 s < t < 1221 s) with a value of approximately 160 K. Minimal temperature differences are seen for t < 300 s and t > 1500 s, corresponding to the completion of the first stage heating ramp and the shutoff cooling phase, respectively.The intermediate temperature differences, corresponding to time sequences forward and aft of the dwell interval, generally show linearly increasing and decreasing trends, respectively.
The evolution of the hydrostatic stress (σ h = 1/3I 1 ) where I 1 is the first stress invariant (I 1 = σ 1 + σ 2 + σ 3 ), at points 1 and 2, is shown in Figure 6(a).The hydrostatic stress was selected over an analysis utilizing the four applicable components of the stress tensor (σ z , σ r , σ θ , and σ rz ) due to its relative simplicity (scalar representation) and, more importantly, its utility in describing changes in volume, which by definition is represented by the first principle stress invariant.This parameter thus helps to account for the processes associated with rates of densification and has been used in like manner in various other powder densification processes [18].As indicated in Figure 6(a), the stress evolutions at points 1 and 2 are nearly mirror images of each other, only differing with respect to magnitude.Further, contrasting Figure 6(a) with Figure 3(a), the characteristic periods of heating, dwell, and cooling are distinguishable.For example, for the dwell period associated with the peak temperature, between t = 857 s and t = 1221 s, the hydrostatic stress for both points 1 and 2 is also shown to be relatively constant, and the final cooling period associated with t > 1516 s shows relatively steep stress gradients back to the original null stress condition.While much of the stress distribution at point 1 is compressive, the dominant influence of the tensile nature of the thermal stresses is evident at time t < 857 s, where the two successive heating ramps are manifested.In contrast, the hydrostatic stress at point 2 is compressive during these early heating stages and tends to have a relatively minimal influence thereafter.
The hydrostatic stress distribution, plotted at time t = 1000 s and as a function of radial distance within the sample, is shown in Figure 6(b).Like the temperature plot of Figure 4(b), the compressive stress is nonhomogeneous and exhibits a significant difference (∼35 Mpa) between the computed values at point 1 and 2. This nonuniform stress distribution, like the results found with respect to temperature, can affect the resulting homogeneity of the microstructure, resulting, for example, in nonuniform grain sizes.

Conclusions
In this work we presented a transient, FAST, finite element, and heat transfer model of aluminum oxide sintering.The model utilized a coupled electrical, thermal, and mechanical analysis approach that included the use of a proportional integral derivative (PID) controller, which adjusted the applied voltage in accordance with the specification of a predetermined heating rate.Comparisons of the model with physical experiments revealed excellent agreement in terms of quantifying the temperature evolution at a fixed point outside of the sample.Inside the sample, the numerical results showed significant nonhomogeneous temperature and hydrostatic stress distributions, which could both potentially compromise the integrity of the material microstructure.
Due to the limiting assumption of constant density within the sample, the results presented here may only be fully applicable at the final stages of sintering, wherein the sample microstructure is fully dense and of negligible porosity.In order to expand the functionality of the present model, future implementations will necessarily include the densification kinetics associated with the changing porosity.This additional coupling will allow for potentially better predictions over the entire sintering cycle.

Figure 1 :
Figure 1: (a) The field assisted sintering technique (FAST) in this study utilizes a tool assembly composed of specialized graphite composite (R7710 [21]).(b) A schematic of the FAST geometry, illustrating the placement of the Al 2 O 3 sample and the centrally located hole(s) for the placement of a thermocouple.

Z
(b).The die is 48 mm tall and has external and internal diameters of 40 mm and 20 mm, respectively.Each of the opposing punches is 35 mm tall and has a diameter of 20 mm.A central hole is drilled in each punch and used for the accommodation of a pyrometer Applied load P(t) Applied V (t) T = 673.15K T = 673.15K

Figure 2 :
Figure 2: (a) The associated thermal, electrical, and mechanical boundary conditions are placed on a simplified 2D axisymmetric geometry.The cylindrical coordinate axis is placed along the central axis with origin at the center of the sample.Points 1, 2, 3, 4, and 5 correspond to various reference locations used throughout this work.The mesh, composed of 12,428 three-node triangular elements, was applied using a Delaunay triangulation.(b) The experimental and simulated applied uniaxial loading as a function of time.

Figure 3 : 2 Figure 4 :
Figure 3: Time evolution profiles of the temperature and the electrical potential taken at point 5, located just above the Al 2 O 3 sample.The temperature profile (a) compares the results from the physical experiments and the simulation.The electrical potential profile (b) shows the corresponding changes in electric potential that resulted from computations via a PID controller and the associated heating/cooling rate inputs.

Figure 3 (Figure 5 :
Figure 5: Temperature contour plots corresponding to four different time periods: t = 300 s, t = 600 s, t = 1000 s, and t = 1300 s.The contours in concert with the isotherms further highlight the strong temperature gradients which exist particularly at times t = 600 s and t = 1000 s.

Figure 6 :
Figure 6: (a) The hydrostatic stress evolution at points 1 and 2 showing the opposing forces due to the uniaxial loading and thermal expansion acting over the specified time intervals.(b) The hydrostatic stress as a function of radial position inside the Al 2 O 3 sample at line 1-2.Similar to the temperature distribution within the sample, the stress is also significantly nonhomogeneous.

Table 1 :
FAST sample and tool dimensions.

Table 2 :
Properties of R7710 (graphite) and alumina as functions of temperature (273.15K