## Abstract

Rapid, energy-efficient sintering of materials comprised heterogeneous powders is of critical importance in emerging technologies where traditional manufacturing processes may be difficult to apply. In particular, electrically aided sintering, which uses the material's inherent resistance to flowing current—resulting in Joule heating to bond the powder components—has great promise because it produces desired materials without much post-processing. Furthermore, it has advantages over other methods, such as high purity of processed materials, in particular, because there are few steps during the approach. In order to electrically process the material properly, one must ascertain the externally applied field to properly Joule heat the various material components in the powder mixture. The Joule-heating field is mathematically expressed by the inner product (** J**⋅

**) of the current (**

*E***) and electric (**

*J***) fields throughout the system. This study develops estimates for the Joule-heating fields carried by each phase in a powder mixture, using knowledge of only the externally applied current, and the material properties of the components comprising the mixture. These estimates are useful in guiding and reducing time-consuming material synthesis involving laboratory experiments and/or large-scale numerical simulation.**

*E*## 1. Introduction

Generally, sintering refers to processing a compacted powder material, which is brittle (‘green’), by directly heating it to 70–90% of the melting temperature, for example, by placing it in a furnace, typically with three chambers: (i) a burn-off chamber to the vaporize lubricants (used for easy pouring and compaction), (ii) a high-temperature chamber to sinter, and (iii) a cooling chamber to ramp down the temperature. The binding occurs by small-scale mechanisms involving diffusion, plastic flow, recrystallization, grain growth and pore shrinkage.^{1} Sintering has distinct advantages over other methods, primarily because of the relatively few number of steps during the overall process (thus retaining the material purity) and the production of a near final shape of the desired product without much post-processing. However, because powder processing is typically more expensive than other material processing involving full-blown melting, research is ongoing to improve the steps in the process. In this regard, electrically aided sintering techniques for heat delivery are promising. The key quantity of interest here is the amount of heat generated from running a current through a material, denoted *H* (a rate), which feeds into the first law of thermodynamics,
1.1
In equation (1.1), *ρ* is the mass density, *w* is the stored energy per unit mass, ** T** is the Cauchy stress,

**is the displacement field,**

*u***is the heat flux and**

*q**H*=

*a*(

**⋅**

*J***) is the rate of electrical energy absorbed owing to Joule heating, where**

*E***is the current,**

*J***is the electric field and 0≤**

*E**a*≤1 is an absorption constant. Our objective in this paper is to determine the phase-wise load-shares of the Joule field, denoted

*H*=

**⋅**

*J***, carried by the components in the heterogeneous powder mixture.**

*E*It is important to realize that heterogeneous mixtures (microstructures) distort the electrical and current field within the material. For electrically aided sintering to be properly controlled, in particular, for heterogeneous powder mixtures, one needs accurate characterizations of the electrical loads carried by each of the phases in the system. In this paper, as a model problem, we will consider a statistically representative volume element (RVE of volume |*Ω*|) of a two-phase dielectric medium, as depicted in figure 1. *We assume that the material has been properly compacted so that there are no gaps between the phases* (*an idealization*). The microscale properties are characterized by a spatially variable electrical conductivity, ** σ**(

**). For such a sample, one can decompose the electrical field carried by each phase in the material as follows: 1.2 the current can be decomposed as 1.3 and the Joule-heating field as 1.4 where is a volume averaging operator;**

*x**v*

_{1}and

*v*

_{2}are the volume fractions of phases 1 and 2, respectively (

*v*

_{1}+

*v*

_{2}=1). We denote

*v*

_{1}(〈

*H*〉

_{Ω1}/〈

*H*〉

_{Ω}) and

*v*

_{2}(〈

*H*〉

_{Ω2}/〈

*H*〉

_{Ω}) as the ‘load-shares’ because 1.5 The objective of this paper is to determine the load-shares as functions of known (

*a priori*) quantities, 1.6 and 1.7 where

*σ*

_{1}and

*σ*

_{2}are the conductivities of phase 1 and phase 2, respectively. The overall volume averages, 〈

**〉**

*E*_{Ω}and 〈

**〉**

*J*_{Ω}, are considered known

*because they can be determined by the boundary values from well-known results (discussed in §2): (i) the Average Electric Field Theorem and (ii) the Average Current Theorem*. The presentation is as follows:

— expressions are developed for the current-field (

) and electric-field (*J*) distribution for each component in the mixture;*E*— expressions are developed for the Joule-heating field current distribution for each component in the materials (

⋅*J*);*E*— bounding principles are used to provide estimates of the overall response of the material;

— asymptotic cases of extreme powder mixtures of insulators and super-conductors are considered;

— simple estimates for the time to heating are provided; and

— extensions, involving numerical methods, are discussed.

### Remark 1.1

The mathematical form for Joule heating can be motivated by taking Faraday's Law
1.8
and Ampere's Law
1.9
where ** D** is the electric field flux,

**is the magnetic field,**

*H***is the magnetic field flux and forming the difference between the inner product of the electric field with Ampere's Law and the inner product of the magnetic field with Faraday's Law, 1.10 where is the electromagnetic energy and**

*B***=**

*S***×**

*E***is the Poynting vector. This relation can be rewritten as 1.11 Equation (1.11) is usually referred to as Poynting's theorem. This can be interpreted,**

*H**for simple material laws, where the previous representation for*, as stating that the rate of change of electromagnetic energy within a volume, plus the energy flowing out through a boundary, is equal to the negative of the total work performed by the fields on the sources and conduction. This work is then converted into thermo-mechanical energy (‘Joule heating’,

*W*holds*H*in equation (1.1)). Joule heating stems from ions being pulled through a medium by electromagnetic fields, which generate heat when they collide with their surroundings.

### Remark 1.2

In this work, we assume that the material has been properly compacted, and we focus solely on the electrical heating aspects of this process. There has been considerable research activity in *non-electrical* compaction of powders, for example, see Brown & Abou-Chedid (1994), Fleck (1995), Tatzel (1996), Akisanya *et al.* (1997), Domas (1997), Anand & Gu (2000), Gu *et al.* (2001), Gethin *et al.* (2003), Zohdi (2003*a*) and Ransing *et al.* (2004).

### Remark 1.3

Generally, for detailed pointwise information, for example, localized effects in the matrix ligaments between particles (‘hot spots’), one needs to solve boundary-value problems posed over a statistically RVE sample of heterogeneous media. This will be discussed at the end of this paper. However, the essential issue is that time-transient effects lead to coupling of electrical and magnetic fields, and the only viable approach is to employ direct numerical techniques to solve for Maxwell's equations. Generally, these equations are strongly coupled. Additionally, if the local material properties are thermally sensitive, and Joule heating is significant, then the first law of thermodynamics must also be solved, simultaneously. Numerical techniques for the solution of coupled boundary-value problems posed over heterogeneous electromagnetic media, undergoing thermo-mechano-chemical effects can be found in Zohdi (2010).

## 2. The controllable quantities: 〈*J*〉_{Ω} and 〈*E*〉_{Ω}

*J*

*E*

For our model problem, two physically important test boundary (∂*Ω*) loading states are notable on a sample of heterogeneous material (figure 1): (i) applied electric fields of the form and (ii) applied current field of the form , where and are constant electric field and current field vectors, respectively. Clearly, for these loading states to be satisfied within a macroscopic body under non-uniform external loading, the sample must be large enough to possess small boundary field fluctuations relative to its size. Therefore, applying (i)- or (ii)-type boundary conditions to a large sample is a way of reproducing approximately what may be occurring in a statistically representative microscopic sample of material within a macroscopic body. The following two results render 〈** J**〉

_{Ω}and 〈

**〉**

*E*_{Ω}as controllable quantities, via the boundary loading.

—

*The Average Electric Field Theorem*. Consider a sample with boundary loading . We make use of the identity 2.1 and substitute this in the definition of the average electric field, 2.2 Thus, if**∇**×=*E***0**, then , when .—

*The Average Current Field Theorem*. Consider a sample with boundary loading . We make use of the identity 2.3 and substitute this in the definition of the average current, 2.4 Thus, if**∇**⋅=0, then when .*J*

### Remark 2.1

The importance of the *Average Electric Field Theorem* and the *Average Current Field Theorem* is that we can consider 〈** E**〉

_{Ω}and 〈

**〉**

*J*_{Ω}to be controllable quantities, via or on the boundary. Applying these boundary conditions should be made with the understanding that these idealizations reproduce what a RVE (which is much smaller than the structural component of intended use) would experience within the system of intended use. Uniform loading is an idealization and would be present within a vanishingly small microstructure relative to a finite-sized engineering (macro)structure. These types of loadings are somewhat standard in computational analyses of samples of heterogeneous materials (see Zohdi 2003

*b*, 2010; Zohdi & Wriggers 2008; Ghosh 2011; Ghosh & Dimiduk 2011).

### Remark 2.2

In the analysis that follows, we will use the following energy–power relation: 2.5 which is referred to as an ergodicity condition in statistical mechanics (Kröner 1972; Torquato 2002) and as a Hill-type condition in the solid mechanics literature (Hill 1952). This is essentially a statement that the microenergy (power) must equal the macroenergy (power). Equation (2.5) is developed by first splitting the current and electric fields into mean (average) and purely fluctuating (zero mean) parts. For the current field, one has , where and for the electric field , where . The product yields 2.6 because and . The ergodicity assumption is that , as the volume (relative to the inherent length scales in the microstructure). The implication is that, as the sample becomes infinitely large, is purely fluctuating and hence . In other words, the product of two purely fluctuating random fields is also purely fluctuating. These results are consistent with the use of the uniform boundary loadings introduced earlier because they can be shown to satisfy equation (2.5).

## 3. Concentration tensors and load-shares

A useful quantity that arises in the analysis of heterogeneous dielectric materials is the effective conductivity, ** σ***, defined via

^{2}3.1 which is the ‘property’ (a relation between volume averages) used in micro–macroscale analyses. Decomposing the left-hand side yields 3.2 where 3.3

*C*_{E,2}is known as the electric field concentration tensor. Thus, the product of

*C*_{E,2}with 〈

**〉**

*E*_{Ω}yields 〈

**〉**

*E*_{Ω2}.

*It is important to realize that once either*

*C*_{E,2}

*or*

*****

*σ**are known, the other can be computed*.

In order to determine the concentration tensor for phase 1, we have from equation (1.2) 3.4 where 3.5 Note that equation (3.5) implies 3.6

Similarly, for the current, we have 3.7 Thus, 3.8 and 3.9 where 3.10 We remark that equation (3.10) implies 3.11 Summarizing, we have the following results:

### Remark 3.1

As a consequence of previous results, the Joule fields can be written in a variety of useful forms, 3.12

## 4. Joule-heating load-shares

Using equation (3.12), the Joule fields can be bounded as follows (using the Cauchy–Schwarz inequality): 4.1

If the overall property is isotropic, and each of the constituents is isotropic (for example, a microstructure comprised of an isotropic binder (fine-scale powder) embedded with randomly distributed isotropic particles), then we have the following, *C*_{E,i}=*C*_{E,i}**1** where, for a two-phase material,
4.2
and *C*_{J,i}=*C*_{J,i}**1**, leading to
4.3
Thus, in the case of isotropy, equation (4.1) asserts
4.4
The products of the concentration functions take the following forms:
4.5
and
4.6
Because the concentration functions depend on *σ**, which in turn depends on *σ*_{1}, *σ*_{2}, *v*_{2} and the microstructure, we need to employ estimates for *σ**. One class of estimates are the Hashin–Shtrikman bounds (Hashin & Shtrikman 1962) for two isotropic materials with an overall isotropic response,
4.7
where the conductivity of phase 2 (with volume fraction *v*_{2}) is larger than phase 1 (*σ*_{2}≥*σ*_{1}). Provided that the volume fractions and constituent conductivities are the only known information about the microstructure, the expressions in equation (4.7) are the tightest bounds for the overall isotropic effective responses for two-phase media, where the constituents are both isotropic. A critical observation is that the lower bound is more accurate when the material is composed of high-conductivity particles that are surrounded by a low-conductivity matrix (denoted case 1) and the upper bound is more accurate for a high-conductivity matrix surrounding low-conductivity particles (denoted case 2).

The previous comments on the accuracy of the lower or upper bounds can be further qualitatively explained by considering the two cases with 50 per cent low-conductivity material and 50 per cent high-conductivity material. A material with a continuous low-conductivity (fine-scale powder) binder (50%) will isolate the high-conductivity particles (50%), and the overall system will not conduct electricity well (this is case 1 and the lower bound is more accurate), while a material formed by a continuous high-conductivity (fine-scale powder) binder (50%) surrounding low-conductivity particles (50%, case 2) will, in an overall sense, conduct electricity better than case 1. Thus, case 2 is more closely approximated by the upper bound and case 1 is closer to the lower bound. Because the true effective property lies between the upper and lower bounds, one can construct the following approximation:
4.8
where 0≤*ϕ*≤1. *ϕ* *is an unknown function of the microstructure.* However, the general trends are (i) for cases where the upper bound is more accurate, , and (ii) for cases when the lower bound is more accurate, . Explicitly, for the product of concentration functions, embedding the effective property estimates, we have
4.9
and
4.10

### Remark 4.1

There are a vast literature of methods, dating back to Maxwell (1867, 1873) and Rayleigh (1892), to estimate the overall macroscopic properties of heterogeneous materials. For an authoritative review of (i) the general theory of random heterogeneous media, see Torquato (2002), (ii) for more mathematical homogenization aspects, see Jikov *et al.* (1994), (iii) for solid mechanics inclined accounts of the subject, see Hashin (1983), Mura (1993) and Nemat-Nasser & Hori (1999), (iv) for analyses of cracked media, see Sevostianov *et al.* (2001), and (v) for computational aspects, see Ghosh (2011), Ghosh & Dimiduk (2011) and Zohdi & Wriggers (2008). Tighter estimates, including generalized N-phase bounds can be found in Torquato (2002).^{3}

### Remark 4.2

The governing equation used in developing effective conductivity bounds is **∇**⋅** J**=0, which stems from taking the divergence of Ampere's Law:

**∇**⋅(

**∇**×

**−∂**

*H***/∂**

*D**t*−

**)=0, one obtains, because**

*J***∇**⋅(

**∇**×

**)=0, 4.11 where is the charge per unit volume. Thus, if , then**

*H***∇**⋅

**=0. If one employs the constitutive relation**

*J***=**

*J***⋅**

*σ***, then this allows for Hashin–Shtrikman type estimates to be used for the effective conductivity, as does**

*E***∇**⋅

**=0 (which is valid only when ) for estimates of the effective permittivity, 〈**

*D***〉**

*D*_{Ω}=

***⋅〈**

*ϵ***〉**

*E*_{Ω}, when

**=**

*D***⋅**

*ϵ***. For example, one case when these two physical situations are compatible is when**

*E***=**

*E*

*σ*^{−1}⋅

**=**

*J*

*ϵ*^{−1}⋅

**⇒**

*D***=(**

*J***⋅**

*σ*

*ϵ*^{−1})⋅

**.**

*D*## 5. Examples of Joule-heating load-sharing

### (a) A general dielectric mixture

Figures 2 and 3 illustrate a surface (using ) in parameter space (*σ*_{1}/*σ*_{2},*v*_{2}) for the normalized Joule-heating load-share, *v*_{i}(〈*H*〉_{Ωi}/〈*H*〉_{Ω}), of each component, *i*=1,2. The plots illustrate the proportion of the Joule heating that will be delivered to each phase in the system. Directly from equations (4.9) and (4.10), the load-share quantities of interest are
5.1
and
5.2
The trends are

— for phase 1: decreasing the volume fraction of phase 2 (

*v*_{2}), for fixed*σ*_{2}/*σ*_{1}, leads to a larger load-share for phase 1, whereas decreasing the mismatch*σ*_{2}/*σ*_{1}, for a fixed*v*_{2}, leads to an increased load-share for phase 1, for a fixed volume fraction*v*_{2}and— for phase 2: increasing the volume fraction of phase 2 (

*v*_{2}), for fixed*σ*_{2}/*σ*_{1}, leads to a larger load-share for phase 2, whereas increasing the mismatch*σ*_{2}/*σ*_{1}, for a fixed*v*_{2}, leads to an extremely slight change in the load-share of phase 2 (it is virtually flat).

### (b) An extreme mixture: high-conductivity (‘superconducting’) particles in a low-conductivity matrix

For the case of high-conductivity particles (phase 2) in a lower conductivity matrix (phase 1), we have
5.3
Inserting this expression into the Hashin–Shtrikman bounds and taking the limit as yields (*σ*_{2} tending to infinity)
5.4
where the lower (Hashin–Shtrikman) bound is more accurate (). Correspondingly, for the concentration tensors for phase 1 (assuming isotropy),^{4}
5.6
and for phase 2 (particle),
5.7
Forming the products yields
5.8
and
5.9
The expressions are appropriate for small *v*_{2} (superconducting particles in a binding matrix). Thus, we have for the load-share
5.10
while for phase 2 (particle superconductor, no Joule field),
5.11

### Remark 5.1

As (no particle material; phase 2), the expressions collapse to restrictions on the pure matrix (here, phase 1) material.

### (c) An extreme mixture: low-conductivity (‘insulator’) particles in a high-conductivity matrix

For the case of low-conductivity particles (phase 1) in a high-conductivity matrix (phase 2), we have
5.12
Inserting this expression into the Hashin–Shtrikman bounds and taking the limit as (*σ*_{1} tending to zero) yield
5.13
where the upper (Hashin–Shtrikman) bound is more accurate (). Correspondingly, for the concentration tensors (as ), for phase 1 (particle),
5.14
and for phase 2 (matrix),
5.15
The expressions are appropriate for large *v*_{2} (insulating particles in a binding matrix),
5.16
and
5.17
Thus, we have for the load-shares, for phase 1 (particle insulator, no Joule field),
5.18
and for phase 2 (matrix),
5.19

### Remark 5.2

As (no particle (here, phase 1) material), the expressions collapse to restrictions on the pure matrix (here, phase 2) material.

## 6. Conclusions and extensions

Short of large-scale numerical simulations, one can make rough estimates for the time scales for heating the components materials to reach a target temperature by *ignoring stress–power and conduction in the first law* (equation (1.1)). If we further assume that the temperature is uniform in each phase with for each component material (*i*), we have
6.1
which can be integrated to
6.2
where Δ*θ*=*θ*(*t*)−*θ*(*t*=0). Specifically, for a two-phase material,
6.3
and
6.4
yield, with *σ**≈*ϕσ*^{*,+}+(1−*ϕ*)*σ*^{*,−}, for phase 1, the time to heat the material to the desired level (〈Δ*θ*〉_{Ω1}),
6.5
and for phase 2 (〈Δ*θ*〉_{Ω2}),
6.6
Clearly, once the design parameter estimates have been made to estimate the processing time, more detailed information, for example, localized effects in the matrix ligaments between particles (‘hot spots’), can be generated via only numerical simulation. To determine the generation of transient electromagnetic fields, temperature fields, stress fields (owing to both Joule heating and electromagnetically induced body forces) and chemical fields, this requires the solution to the time-transient forms of (i) Maxwell's equations, (ii) the first law of thermodynamics, (iii) the balance of linear momentum, and (iv) reaction–diffusion laws. In order to accurately capture the coupled (transient) electromagnetic, thermal, mechanical and chemical behaviour of a complex material, Zohdi (2010) addressed the modelling and simulation of such strongly coupled systems using a staggered temporally adaptive finite difference time domain (FDTD) method. Of particular interest was to provide a straightforward modular approach to finding the effective dielectric (electromagnetic) response of a material, incorporating thermal effects, arising from Joule heating, which alter the pointwise dielectric properties such as the electric permittivity, magnetic permeability and electric conductivity. Because multiple field coupling is present, a staggered, temporally adaptive scheme was developed to resolve the internal microstructural electric, magnetic and thermal fields, accounting for the simultaneous pointwise changes in the material properties. This approach also incorporated the coupled chemical and mechanical fields. We remark that there are a variety of computational electromagnetic methods. The most widely used technique is the FDTD, which is ideally suited to the problems of interest in this work. However, there are other methods, such as (i) *the Multi Resolution Time Domain Method*, which is based on wavelet-based discretization, (ii) *the Finite-Element Method*, which is based on discretization of variational formulations and which is ideal for irregular geometries,^{5} (iii) *the Pseudo Spectral Time Domain Method*, which is based on Fourier and Chebyshev transforms, followed by a lattice or grid discretization of the transformed domain, (iv) *the Discrete Dipole Approximation*, which is based on an array of dipoles solved iteratively with the conjugate gradient method and a fast Fourier transform to multiply matrices, and (v) *the Method of Moments*, which is based on integral formulations employing boundary-element method discretization, often accompanied by the fast multi-pole method to accelerate summations needed during the calculations, and (vi) *the Partial Element Equivalent Circuit Method*, which is based on integral equations that are interpreted as circuits in discretization cells.

In Zohdi (2010), FDTD was combined with a staggering solution framework to solve the coupled dielectric material systems of interest. The general methodology is as follows (at a given time increment): (i) each field equation is solved individually, ‘freezing’ the other (coupled) fields in the system, allowing only the primary field to be active and (ii) after the solution of each field equation, the primary field variable is updated, and the next field equation is treated in a similar manner. For an *implicit* type of staggering, the process can be repeated in an iterative manner, while for an *explicit* type, one moves to the next time step after one ‘pass’ through the system. As the physics changes, the field that is most sensitive (exhibits the largest amount of relative non-dimensional change) dictates the time-step size. Because the internal system solvers within the staggering scheme are also iterative and use the previously converged solution as their starting value to solve the system of equations, a field that is relatively insensitive at given stage of the simulation will converge in very few internal iterations (perhaps even one). The overall goal was to deliver solutions where the staggering (incomplete coupling) error is controlled and the temporal discretization accuracy dictates the upper limits on the time-step size. Generally speaking, the staggering error, which is a function of the time-step size, is time dependent and can become stronger, weaker or possibly oscillatory, and is extremely difficult to ascertain *a priori* as a function of the time-step size. Therefore, to circumvent this problem, an adaptive staggering strategy was developed to provide accurate solutions by iteratively adjusting the time steps. Specifically, a sufficient condition for the convergence of the presented fixed-point scheme was that the spectral radius (contraction constant of the coupled operator), which depends on the time-step size, must be less than unity. This observation was used to adaptively control the time-step sizes while simultaneously controlling the coupled operator's spectral radius, in order to deliver solutions below an error tolerance within a prespecified number of desired iterations. This recursive staggering error control can allow for substantial reduction of computational effort by the adaptive use of large time steps, when possible. Furthermore, such a recursive process has a reduced sensitivity (relative to an explicit staggering approach) to the order in which the individual equations are solved because it is self-correcting. For more details, see Zohdi (2010). The further development of numerical methods for electrically aided sintering simulation is under further investigation by the author.

## Acknowledgements

The author expresses his gratitude to Ms Cora Schillig and Mr Gary Merrill of the Siemens corporation for support of this research.

## Footnotes

↵1 An oxygen-free environment is best to minimize oxides.

↵2 Implicitly, we assume that (i) the contact between the phases is perfect and (ii) the ergodicity hypothesis is satisfied (see Kröner (1972) or Torquato (2002)).

↵3 Such N-phase bounds go well beyond the simple Wiener bounds (Wiener 1910), .

↵4 These expressions are asymptotically consistent with the identities 5.5

↵5 In particular, see Demkowicz (2006) and Demkowicz

*et al.*(2007) for the state-of-the-art in adaptive finite-element methods for Maxwell's equations.

- Received December 31, 2011.
- Accepted February 8, 2012.

- This journal is © 2012 The Royal Society