# BAKING THE CARBON CAKE Homogenisation Applied to the Electrical Calcination of Carbon Materials

By Caoimhe Rooney

The metallurgy industry is a gold-mine for mathematical treasures, from partial differential equation models describing crust formation (which notoriously hinders furnace operation), to population balance equations explaining the formation and growth of microsilica particles (which impacts the performance of composite materials). Through the collaboration of the University of Oxford and Norwegian companies NORCE, Elkem and Eramet, mathematicians strive to tackle such industrial challenges by developing and solving mathematical models. A common feature of many metallurgical processes is the use of electrical current flow to heat a granular material with the desire to manipulate its chemical and physical properties. Understanding and optimising this phenomenon forms the basis of my PhD research.

Calcination describes the heat treatment of carbon in a furnace to produce a partially-graphitised material for use in electrodes and for other metallurgical applications. The carbon is in the form of anthracite particles which are irregular in size and rather angular in shape, and make only a few microscopic contacts with their neighbours. Figure 1 shows a CT scan of such particles. An electric current is passed through the ‘coke bed’ of anthracite particles, causing Ohmic heating and high temperatures which result in the chemical and structural transformation of the material. To explain behaviour in these furnaces, one must understand how electric current and thermal energy propagate through the granular material within it. This is a complex process due to the irregular and inhomogeneous shapes of the particles, with a particular dependence on the nature of the contacts between particles. These contacts are microscopic in comparison with the scale of the furnace, but their impact on the electrical and thermal conductivity is critical. Therefore, understanding the process on a microscopic level is imperative to understanding the macroscopic behaviour of the furnace. Figure 1: CT scan slices of cylindrical crucible containing carbon particles. The left image shows the axial slice corresponding to the vertical position indicated in the right image

Through experimental exploration and scaling arguments we have concluded that the electrical resistance at the contacts between particles contributes substantially to the bulk resistivity of the coke bed. Similarly, a sound understanding of thermal transfer between and within particles is imperative to fully explain the temperature distribution of the material. Understanding the behaviour of such mechanisms on the scale of a single particle is often dealt with through the use of computational models such as DEM (Discrete Element Methods). However, because of the great discrepancy between the length scale of the particles and the length scale of the furnace, we can exploit asymptotic homogenisation theory to simplify models.

To illustrate the efficacy of mathematical homogenisation, consider steady electrical current conduction within a simplified representation of the coke bed of unit height, where we model particles as regular squares of dimension $\delta$, $\ll$ 1 as in Figure 2. We capture the effect of gasesous voids by enforcing insulating boundary conditions on a portion of the particle boundary. The current flow induces Ohmic heating within the material, therefore we model electrical and thermal conduction by \begin{aligned} \nabla \cdot(\sigma(T)\nabla\phi(\mathbf{x})) &=0, &&\quad \mathbf{x}\in \Omega, \\ \nabla \cdot(k(T)\nabla T(\mathbf{x})) + \sigma(T) \nabla\phi(\mathbf{x})^\top\nabla\phi(\mathbf{x}) &=0, &&\quad \mathbf{x}\in\Omega, \end{aligned}

where $\phi$ is the electric potential, $T$ is the temperature, $x$ is the spatial variable and the electrical and thermal conductivities, which depend on the local temperature, are denoted $\sigma(T)$ and $k(T)$ respectively. We must therefore solve a coupled, non-linear system of partial differential equations.

The domain is subject to mixed Dirichlet-Neumann exterior boundary conditions for both $\phi$ and $T$. Internally, we enforce two types of simplified boundary conditions upon each particle to differentiate between particle-to-particle and particle-to-gas interfaces. Assume that current and heat conduct continuously at contact interfaces, hence we enforce a continuity condition where the particles make contact with each other. Figure 2: A simplified two-dimensional representation of a collection of square particles. Each particle contacts its neighbours at an interface of width $\epsilon$, and the remaining particle boundary is assumed to be insulating. We enforce mixed Dirichlet-Neumann boundary conditions on the exterior of the domain.

These contacts are approximated as a gap of width $\epsilon$ $\ll$ $\delta$ in the boundary between particles. For simplicity, we neglect thermal radiation and assume that current cannot conduct through gas, therefore the remaining interface between the particles is assumed to be an insulating boundary. This indicates that electrical and thermal transfer is made via particle-to-particle contacts only, and not through the gas. Hence the boundary conditions are \begin{aligned} \phi(\mathbf{x}) = T(\mathbf{x}) &= 1, &&\quad \mathbf{x}\in\partial\Omega_1,\\ \phi(\mathbf{x}) = T(\mathbf{x}) &= 0, &&\quad \mathbf{x}\in\partial \Omega_2,\\ \nabla\phi(\mathbf{x})\cdot\mathbf{n} = \nabla T(\mathbf{x})\cdot\mathbf{n} &= 0, &&\quad \mathbf{x}\in\partial\Omega_3, \end{aligned}

where $\partial\Omega_i, \, i = 1,2,3$ represent the boundaries of the domain, and more specifically, $\partial\Omega_3$ includes the insulating boundaries between the particles.

To solve this problem numerically, one must use a very fine mesh to capture the behaviour of the solution within each microscopic particle. Moreover, the mesh must be fine enough to describe the behaviour of the potential near to the contact interfaces, since this is crucial when simulating current conduction through granular material. Given that the dimension of the particle contact is orders of magnitude smaller than the bulk dimension $\epsilon$ $\ll$ $\delta$ $\ll$ l, this constitutes an extremely fine mesh. The problem is therefore extremely computationally expensive, and as we increase the complexity of the problem by introducing more realistic features such as a third dimension and different particle geometries, a numerical solution obtained in this way becomes unfeasible.

Due to the repetitive or periodic nature of the microstructure of the domain, and
given that $\delta$ $\ll$ l, we can simplify the problem using the method of homogenisation. Put simply, this method allows us to smooth out the inhomogeneity of the microstructure by capturing its impact on the electric potential through an effective conductivity in a simplified macroscopic problem, defined on a homogeneous domain, which can be solved using a significantly coarser mesh than the full problem. This effective-medium problem captures the important features of the problem whilst minimising the complexity.

(a)                                    (b)                                   (c)

(a) Numerical Solutions solved for $\delta$ = 0.1 $\epsilon$ = 0.004
(b) Numerical Solutions solved for $\delta$ = 0.05 $\epsilon$ = 0.002
(c) Homogenised Solution

Figure 3: Numerical solution 1-5 of the electric potential, solved directly and by homogenisation. $\delta$ $\delta$ $0.1$ $0.1$ $0.005$ $0.005$ $0.025$ $0.025$ Relative L2 error $0.037$ $0.037$ $0.022$ $0.022$ $0.013$ $0.013$ Computational time(s) $14.28$ $14.28$ $59.91$ $59.91$ $274.8$ $274.8$ $\sim 400 \%$ $\sim 400 \%$ $\sim 1650\%$ $\sim 1650\%$ $\sim 7590 \%$ $\sim 7590 \%$

Table 1: Summary of error between numerical solution and homogenised solution for different values of $\delta$, along with the computational time of the numerical solution as a percentage of the homogenisation time ( $3.619$ $s$)

We can solve the problem (1)-(5) for different values of $\delta$ and $\epsilon$ to compare with the homogenised solution. The heatmap of electric potential $\phi$ for each simulation is given in Figure 3. Note that this system requires the electrical and thermal conductivities as inputs. For simplicity, we allow these to be the same, $\sigma(T) =$ and $k(T)$ To capture a trend similar to that observed experimentally for the electrical conductivity of carbon, in this simulation we use $\sigma(T) = 0.05 + \frac{\mathrm{exp}(5(3T-1))}{2(\mathrm{exp}(5(3T-1))+1)}.$

The homogenised solution is a valid approximation of the numerical solution for $\delta$ $\rightarrow$ 0, as indicated by the $\L^2$ error values in Table 1. We note that the homogenisation solution was obtained 4 times faster than a numerical solution of 50 particles, and nearly 76 times faster than a numerical solution of 800 particles. In reality, we wish to model systems containing millions of particles, and in such cases the computational advantage of homogenisation becomes apparent.

Our asymptotic and experimental analysis have indicated that the microscopic contacts contribute significantly to the bulk resistivity of the coke bed. In the example explained here, we have included the effect of constriction resistance, that is, the resistance incurred by current passing through extremely narrow gaps. This assumes a perfect, continuous contact area, however, in reality, we expect surface asperities and oxide layers to interfere with current conduction. We can model such imperfect interfaces with jumps in the potential at contacts between particles. This increases the numerical complexity, as we must consider discontinuous functions.

Furthermore, to include non-regular configurations of particles, we can model the particles as an irregular network of resistors, whose values capture contact and material resistances. By considering Kirchhoff’s and Ohm’s laws, and through the use of the lattice Laplacian and Green’s function, we hope to solve the resistor network problem to gain insight into how randomness and irregularity affects the bulk resistivity of a bed of carbon particles.

Acknowledgments ~ This work is supervised by Colin Please and Sam Howison of the University of Oxford, supported by the EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1), and is part of the research project Electrical Conditions and their Process Interactions in High Temperature Metallurgical Reactors (ElMet/47791). ElMet is a cooperation between NORCE (Project Owner), The Norwegian University of Science and Technology (NTNU) and the industrial companies Elkem and Eramet. The University of Santiago de Compostela, University of Oxford are research partners. The project is financially supported by The Research Council of Norway and the two companies.