A Simplified Bubble Nucleation, Growth and Coalescence Model for Coke Production Process

Abstract

Coke production process, which involves gas generation, cross-linking, gas foaming as well as solidification phenomena, is extremely complex and is difficult to model them without simplification. In this study, a simple phenomenological model was developed based on a gas foaming simulation model of polymer foaming, which enabled us to simulate the number density of bubbles and their distributed size. The model was extended to account their kinetics of bubble nucleation, growth and coalescence in non-isothermal chemical reactions of gas generation and cross-linking. The bubbles size and morphology obtained from the numerical simulation of two different coals agreed with the pictures of the experiment qualitatively. Five different coals were investigated to understand the relationship between the kinetic and final morphology of coke.

1. Introduction

Coke is one of the most important raw materials of pig iron, which was produced by reduction of iron ore. The transformation of coal to coke involves chemical reaction, mass transfer and phase separation, which produces many bubbles in coke.^1,2,3,4) Figure 1 shows the chemical and transport rate processes contributing to volatiles evolution during pyrolysis of softened coal particles. The pyrolysis of coal produces gas, metaplast and char. The produced gas is first dissolved in metaplast. When the gas concentration in the metaplast reached certain level, bubble nucleation occurs to form bubbles. Then gas is transported to bubbles or passed to atmosphere.

Fig. 1.

Chemical and transport rate processes contributing to volatiles evolution during pyrolysis of a softened coal particle.⁴⁾

The size and number of bubbles in other words, porosity of coke determines the important specification of coke which are mechanical strength and permeability of gas through coke. To control the porosity of coke, it is essential to understand how coal is changed in the coke production process. However, direct observation of coke production process is hardly possible because of too high temperature and dangerous atmosphere. Numerical simulation of coke production process is one of the ways to understand the complex phenomena in the process. Moreover, the simulation of coke to draw a final porous structure could be used for further numerical simulation of gas permeation or mechanical characteristics. The numerical simulation may reduce the trial-and-error for optimization of coal formulation among several kinds of coals to obtain coke for production.

The coke production process involves bubble nucleation, growth and coalescence as mentioned above. In the field of polymer foaming community, the bubble nucleation and growth phenomena in polymer matrix are also important to produce polymer foam products. Thus the studies of bubble nucleation and growth with the aid of numerical simulation as well as direct observation have been investigated continuously.^{5,6,7,8,9,10,11)} The industrial polymer foaming process is divided into two categories where chemical and physical foaming.

Chemical foaming is induced by pyrolysis of chemical blowing agent in polymer matrix. Chemical blowing agent produces gas which dissolves in matrix. When the gas concentration reached a certain level, the bubble nucleation starts in polymer matrix.

Physical foaming requires gas dissolution at high gas pressure first. When the pressure is released to atmospheric level, the dissolved gas molecules form nucleus of bubble, which is super saturated state. One of the nuclei is expanded by two driving forces: (1) pressure difference between bubble pressure and atmosphere, and (2) mass transfer of dissolved gas into the bubble.

For both chemical and physical foaming, when two bubbles approach together, they start to coalesce. The bubble coalescence time depends on the viscosity of matrix.⁸⁾ The higher matrix viscosity is, the longer the bubble coalescence time is. The driving force of bubble coalescence is the pressure difference between two bubbles. If the mass transfer is neglected between two bubbles, the pressure difference squeezes polymer matrix from one side to another side, makes it thinner and eventually, ruptures the wall of bubbles.

The gas foaming phenomena of bubble nucleation, growth and coalescence in polymer foaming is widely studied employing numerical simulation and visual observation. Although study of chemical foaming and effect of cross-linking on foaming is limited, the current studies of physical foaming could contribute to understand and develop a kinetic model of coke production process.

In this study, the kinetic model of polymer foaming is extended to the coke production process which involves non-isothermal chemical reaction and bubble coalescence. The bubble growth model is extended to account the change of temperature on the bubble pressure. A simple rule-based model is proposed for complex phenomena of bubble coalescence. A set of parameters of molten state coal is input to the developed model to simulate the bubble size during the coke production process. The result indicated that the snapshot of the model corresponded to a picture of coke in experiment. In the rest of this paper, the theory of model is explained, results of numerical simulation are shown and discussed. Finally, conclusion and remarks are expressed.

2. Theory

The developed model is based on the bubble nucleation and growth model of polymer foaming proposed by Taki.⁷⁾ His model consists of bubble nucleation derived from the classical nucleation theory and bubble growth based on the integration of mass transfer and fluid dynamics. Toramaru used similar expression of bubble nucleation and growth in the expansion process of volcano.¹²⁾ In this study, Taki’s model was extended to non-isothermal condition and gas generation system.

2.1. Bubble Nucleation

The classical nucleation theory is based on water drop formation in supersaturated steam and is widely applied to several nucleation phenomena.¹³⁾   

$$J=f_0{\left[\frac{2\gamma }{\pi m}\right]}^{1/2}exp\left[-\frac{16\pi {\gamma }^{3}F}{3k_{\mathrm{B}}T\left(\overline{c}/k_{\mathrm{H}}-P_{\mathrm{C}}^{*}\right)}\right]\overline{c}{N}_{\mathrm{A}}$$(1)

where J is bubble nucleation rate, γ is surface tension of polymer, π is the circular constant, m is mass of gas molecule, k_B is Boltzmann constant, T is absolute temperature, c is average gas concentration dissolved in matrix, k_H is solubility parameter of gas, $P_{\mathrm{C}}^{*}$ is pressure outside of nucleus gas bubble, N_A is Avogadro number, f₀ is correction factor of probability of bubble nucleation, F is correction factor of free energy barrier of bubble nucleation. When the bubble nucleation rate becomes larger than certain level, it is considered that a number of bubbles are newly formed and placed randomly in the simulation box without overlapping each other.

2.2. Bubble Growth

The diffusion-induced bubble growth model was studied by several integral models of coupling of mass transfer and momentum equations.^14,15,16,17) In this study, the bubble growth model was based on the Rosner and Epstein integration model.¹⁴⁾ They coupled the force balance equation of pressure difference between gas pressure and pressure of matrix and mass transfer of single bubble. Equation (2) is the force balance equation around the bubble.   

$$\frac{dR}{dt}=\frac{R}{4\eta }\left(P_{\mathrm{D}}-P_{\mathrm{C}}-\frac{2\gamma }{R}\right)$$(2)

where R is radius of bubble, t is time, η is viscosity, P_D is bubble pressure, P_C is pressure outside of gas bubble.

The driving force of bubble growth or bubble expansion is the pressure difference between the equilibrium pressure, $P_{\mathrm{D}}^{*}$ calculated by the Laplace equation described by Eq. (3) and bubble pressure, P_D. As describing in Eq. (4), if P_D is larger than $P_{\mathrm{D}}^{*}$ , the bubble is expanded. If P_D is smaller than $P_{\mathrm{D}}^{*}$ , the bubble is shrunk.   

$$P_{\mathrm{D}}^{*}=P_{\mathrm{C}}+\frac{2\gamma }{R}$$(3)

  

$$\frac{dR}{dt}=\frac{R}{4\eta }\left(P_{\mathrm{D}}-P_{\mathrm{D}}^{*}\right)$$(4)

The bubble pressure was determined by the Rosner-Epstein model. The bubble pressure kinetics was derived from the time derivative of equation of state of gas:   

$$\frac{d{n}_g}{dt}=\frac{d}{dt}\left(P_{\mathrm{D}}\frac{4\pi R^3}{3}\frac{1}{\Re T}\right)={4\pi R^{2}D\frac{\partial c}{\partial r}|}_{r=R}$$(5)

where n_g is a number of gas molecules in molar unit, $\Re$ is gas constant, D is diffusion coefficient, c concentration of gas dissolved in matrix.

The rhs of Eq. (5) is the flux of gas into a bubble. To calculate the concentration gradient in the term, the diffusion equation of gas (Eq. (6)) in the matrix is needed to be solved.   

$$\frac{\partial c}{\partial t}=D\left(\frac{2}{r}\frac{\partial c}{\partial r}+\frac{{\partial }^{2}c}{\partial r^2}\right)$$(6)

Although the numerical integration of Eq. (6) to calculate the concentration gradient is feasible, the calculation cost is too high due to the steep gradient at bubble’s surface. To avoid the increase of calculation cost, it is reasonable to introduce integration method proposed by Rosner and Epstein. They assumed a concentration profile of following equation:   

$$\frac{c_{\infty }-c}{c_{\infty }-c_{\mathrm{R}}}={\left(1-\frac{r-R}{\delta }\right)}^2$$(7)

where c_∞ is gas concentration in infinity, c_R is gas concentration in the matrix at the surface of bubble, δ is the width of concentration profile. In the coke production process, the main component of gas is methane. The gas concentration is decreased by the diffusion of gas into bubbles while it is increased by the pyrolysis of coal. The unevenness of gas concentration in matrix (coal or coke) is hard to be calculated. In this model, the average concentration, c is used to calculate the concentration gradient at the interface of bubble.   

$$\frac{\overline{c}-c}{\overline{c}-c_{\mathrm{R}}}={\left(1-\frac{r-R}{\delta }\right)}^2$$(8)

Then the equations of (5), (6) and (8) were integrated and expanded to know the change of radius and temperature on the bubble pressure with an approximation. Detailed derivation is found in the reference.¹⁴⁾   

$$\begin{array}{l}\frac{d{P}_{\mathrm{D}}}{dt}=sgn\left(\overline{c}-c_{\mathrm{R}}\right)\frac{6D{\left(\Re T\right)}^2{\left(\overline{c}-c_{\mathrm{R}}\right)}^{2}R}{P_{\mathrm{D}}{R}^3/T-P_{\mathrm{D}0}{R}_0^3/T_0}\\ \hspace{1em}\hspace{1em}-\frac{3P_{\mathrm{D}}}{R}\frac{dR}{dt}+\frac{P_{\mathrm{D}}}{T}\frac{dT}{dt}\end{array}$$(9)

where P_D0 and R₀ are the bubble pressure and radius of nucleus. Notice that the Eq. (7) is parabolic function. The sgn (c–c_R) gives “+” or “–” depending on the difference of c and c_R in Eq. (9). The first term of rhs describes the mass transfer. The mass transfer totally depends on the diffusion coefficient and c–c_R. The second and third term describe the change of radius and temperature on the bubble pressure, respectively.

2.3. Bubble Coalescence

The direct numerical simulation of bubble coalescence employing explicit or implicit moving boundary system including Eulerian or Lagrangian schemes is hard to perform for many body and long term simulation. A rule-based algorithm was developed to account the large number of bubble coalescences. The bubble coalescence of two bubbles is a biaxial elongation mode of deformation.⁸⁾ The geometry of two bubbles coalescence is shown in Fig. 2. The interface between bubbles ±h is squeezed by force F along the perpendicular plane of bubble-bubble center lines. The impedance of bubble coalescence is viscosity of matrix. If the viscosity is sufficiently low, the surface tension of matrix is prominent or Capillary number. In coke production process, the viscosity is a key factor to control bubble coalescence due to high-viscosity of matrix.

Fig. 2.

Geometry of two bubbles coalescence.⁸⁾

The viscosity at the biaxial elongation was 6η.¹⁸⁾ The driving force of bubble coalescence was the pressure difference between two bubbles. The convenient time scale of bubble coalescence is described by following equation:   

$$\tau =\frac{6\eta \alpha }{\left|P_{\mathrm{D}1}-P_{\mathrm{D}2}\right|}$$(10)

where P_D1 and P_D2 are the bubble pressures of two bubbles. α is non-linear parameter of viscosity i.e., effect of strain-hardening as the cross-linking structure is involved in the coke.

When two bubbles satisfied the Eq. (11), the bubble coalescence time, Δt_c is started to accumulate. After the accumulated time, Δt_c, exceeds the characteristic time, τ, two bubbles are thought to be coalesced and form one bubble.   

$$\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert <R_1+R_2$$(11)

where x_i and R_i are the position and radius of bubble i. The new position of coalesced bubble is assumed to be as follows:   

$$\mathbf{x}=\left({\mathbf{x}}_1+{\mathbf{x}}_2\right)/2$$(12)

The new bubble pressure was assumed to follow the harmonic average of two bubbles’ pressure:   

$$\frac{1}{P_{\mathrm{D}}}=\frac{1}{2}\left(\frac{1}{P_{\mathrm{D}1}}+\frac{1}{P_{\mathrm{D}2}}\right)$$(13)

The assumptions of position and pressure are needed to be assured by ideal experiment or precise numerical simulation. The new radius was calculated from the mass balance of gas in the bubbles.   

$$\frac{8\pi }{3}{R}^3{\left(\frac{1}{P_{\mathrm{D}1}}+\frac{1}{P_{\mathrm{D}2}}\right)}^{-1}=\frac{4\pi }{3}{R}_1^3{P}_{\mathrm{D}1}+\frac{4\pi }{3}{R}_2^3{P}_{\mathrm{D}2}$$(14)

2.4. Average Concentration

The average concentration was calculated by the total mass balance of gas expressed as follows:   

$$\overline{c}=\frac{c_0+{\int }_0^{t}y(t^{\prime })d{t}^{\prime }-\sum _{i=1}^N\frac{4\pi R_i^3{P}_{\mathrm{D},i}}{3\Re T}}{v_0}$$(15)

where c₀ is initial concentration of gas in matrix, y(t) is gas generation rate, v₀ is volume of matrix. The volume of matrix assumed to be constant. The gas generation rate, diffusion coefficient, viscosity were changed with temperature based on the previous studies.

2.5. Gas Generation Rate

The gas generation occurs when the coal becomes molten state at certain level of temperature. The modeling of coal pyrolysis is studied by several authors.¹⁹⁾ There is a peak of gas generation rate because the limited amount of source of gas generation with increase of temperature. In this study, the Gaussian-type equation was employed to simulate the gas generation rate with the increase of temperature.   

$$y=Aexp\left(-\frac{1}{2}{\left(\frac{T-T_{\mathrm{C}}}{w}\right)}^2\right)$$(16)

where y is gas generation rate, A is a constant of gas generation rate, T_C is peak temperature and w is temperature width of half of peak.

2.6. Viscosity

In the coke production process, the coal is melted when the temperature reached to certain temperature. Then the temperature is increased further as the cross-linking reaction of coal occurs with the increase of temperature. The measurements of viscosity in elevating temperature was studied by Hayashi using a needle penetrometry.²⁰⁾ In this model, there is a peak temperature which gives the minimum viscosity, which is assumed to correspond to the peak temperature of gas generation rate. Equation (17) is the model which shows the change of viscosity assuming the Gaussian type function and the maximum and minimum viscosities are η_max and η_min, respectively.   

$$\frac{{log}_{10}{\eta }_{max}-{log}_{10}\eta }{{log}_{10}{\eta }_{max}-{log}_{10}{\eta }_{min}}=exp\left(-\frac{1}{2}{\left(\frac{T-T_{\mathrm{C}}}{w}\right)}^2\right)$$(17)

Equation (17) could give the change of viscosity according to the change of temperature.   

$$\begin{array}{l}{log}_{10}\eta =\\ {log}_{10}{\eta }_{max}-\left({log}_{10}{\eta }_{max}-{log}_{10}{\eta }_{min}\right)exp\left(-\frac{1}{2}{\left(\frac{T-T_{\mathrm{C}}}{w}\right)}^2\right)\end{array}$$(18)

2.7. Diffusion Coefficient

The diffusion coefficient was obtained from the viscosity because the diffusion coefficient and viscosity are explained by the free volume theory. Based on the Doolittle free volume theory for viscosity^21,22) and Vrentas-Duda theory for diffusion coefficient, mathematically, they are interchanged through the free volume fraction.^23,24,25) Although the Vrentas-Duda theory only provides self-diffusion coefficient, the mutual diffusion coefficient in Eq. (9) is mutual diffusion coefficient. It should include the effect of concentration gradient. In this model, the mutual diffusion coefficient is assumed to be identical value of self-diffusion coefficient. The diffusion coefficient is increased with the increase of temperature up to certain level. Then the temperature passed the peak temperature of viscosity. The cross-linking reaction in the coal decreases the diffusivity of gas in the coal. The diffusion coefficient of small molecules in polymer melt is in order of 10^–9 to 10^–11 m²/s.²⁶⁾ The representative diffusion coefficient was 10^–10 m²/s by taking the middle of 10^–9 and 10^–11 m²/s. The minimum viscosity of the molten coal was set 10⁴ Pa·s. Thus, the conversion factor of viscosity to diffusion coefficient was set to 10^–6 m²/s. The diffusion coefficient was calculated by Eq. (19) with the change of viscosity.   

$$D={10}^{-log\eta }{D}_{\mathrm{c}}$$(19)

2.8. Parameters

To understand the relationship between final morphology and the characteristics of different coals, five different sets of parameters were assigned based on literature²⁰⁾ and experiences of engineers in coke production process. The initial and final viscosity of coal and coke is more than 10⁹ Pa·s.²⁰⁾ The viscosity of 10⁹ Pa·s was too high to integrate the bubble growth model of Eqs. (2) and (9) because the bubble pressure was extremely increased by the steep concentration gradient and small rhs of Eq. (2). The bubble growth calculation was failed by the too large bubble pressure with too small bubble radius. In order to relax this situation, the initial viscosity had to be set much lower viscosity, 10⁶ Pa·s as η_max. The onset of bubble nucleation and the viscosity therein are still not clear in experiment. The actual viscosity maybe in the order of 10⁶ Pa·s, which was reasonable value for bubble growth. The final viscosity was set to the same value of the initial viscosity for the sake of simplicity.

The parameters of Eqs. (16) and (18) of four different coals are shown in Table 1. The parameters of Goonyella was determined based on the data.²⁰⁾ The parameters of (1) K9, (2) High MF, (3) White Haven and (4) ASP were totally guessed to realize the characteristics of model coals. Especially, the gas generation rate for bubble nucleation and growth was hard to estimate. The well-known thermo gravimetric analysis of coal measures the gas generation rate from coal particles. The measured gas generation rate was not for the bubble nucleation and growth. It was diffused out from the particles. In the developed model, the escape gas from the coal is not considered yet.

Table 1. Parameters of coals for calculations of gas generation rate and viscosity.

Coal type	Minimum viscosity log ηmin [Pa·s]	Gas generation rate A×107 [N-mL/s]	Peak temperature TC [°C]	Half of peak width w [K]
(1) K9	6.000	4.50	470	20
(2) Goonyella	4.800	4.90	450	29
(3) High MF	4.699	5.00	440	30
(4) White Haven	4.900	5.59	430	27
(5) ASP	4.000	5.88	375	15

Figure 3 shows the change of viscosity, gas generation rate and diffusion coefficient calculated by Eqs. (16), (18) and (19). The (1) K9 is not changed viscosity and diffusion coefficient while the gas generation rate has a peak at 470°C. (1) K9 is rigid coal (kept hard at elevated temperature) and produces relatively low-gas generation rate. (2) Goonyella, (3) High MF and (4) White Haven shows Gaussian type profiles of viscosity, gas generation rate and diffusion coefficient. Those three coals melt with increase of temperature and become “soft” or low-viscosity at the peak temperature of gas generation rate. Moreover, they produce more gas than (1) K9. The viscosity of (5) ASP starts at 10⁴ Pa·s and increases with temperature while the gas generation rate and diffusion coefficient decrease. (5) ASP is soft or low-viscosity at elevated temperature and produces more gas at low temperature.

Fig. 3.

Changes of viscosity, gas generation rate and diffusion coefficient calculated by Eqs. (16), (18) and (19). The vertical straight lines and open circles on viscosity curves indicate the onset of bubble nucleation calculated by the numerical simulation. (Online version in color.)

2.9. Numerical Simulation

The numerical integration of above equations was performed using the Runge-Kutta scheme of MATLAB^® function (ode45). The step-by-step calculation is explained in the flowchart shown in Fig. 4. (1) The temperature, gas generation rate, viscosity and diffusion coefficient were calculated. (2) The bubble nucleation rate was calculated. If the nucleation rate is greater than 1 in Δt_N and in L₀³, it is considered that the bubble nucleation occurs in the model. The number of bubbles is placed in the simulation box randomly without overlapping each other. (3) The bubble growth equations were integrated for each bubble to calculate the radius and pressure of bubble. (4) the distance between every pair of bubbles are calculated. If the distance is shorter than the sum of radius, the bubble coalescence time is accumulated. When the bubble coalescence time becomes larger than the characteristic time of bubble coalescence, the two bubbles were considered to be coalesced. The new bubble pressure and radius are calculated by Eqs. (13) and (14). Table 3 shows a list of parameter to perform the numerical simulation.

Fig. 4.

Flowchart of numerical simulation.

Table 3. Parameters of numerical simulation.

Description	Symbol	Value	Unit
Time scale	tbase	1×10–6	s
Pressure scale	Pbase	101.3×103	Pa
Length scale	Rbase	1×10–9	m
Temperature scale	Tbase	300	K
Gas constant	$\Re$	8.314	J/(mol·K)
Bubble nucleation parameter	f0	1×10–38	–
Mass of gas molecule (methane)	m	2.66×10–23	kg
Wetting factor of bubble nucleation	F	0.1	–
Boltzmann constant	kB	1.38×10–23	J/K
Non-linearity parameter	α	0.1	–
Maximum viscosity	ηmax	106	Pa·s
Conversion factor from viscosity to diffusion coefficient	Dc	10–6	m2/s
Initial size of simulation cell	L03	1×1×1×10–9	m3
Time step of bubble growth	ΔtG	tbase	s
Time step of bubble nucleation	ΔtN	1	s
Initial temperature	T0	650	K
Heating rate	dT/dt	3	K/min
Simulation time	tend	57	min
Surface tension	γ	10×10–3	N/m
Solubility parameter	kH	1×10–7	mol/m3/Pa
Atmospheric pressure	PC	101.3×103	Pa

The result of numerical simulation is a set of positions and size of bubbles. To visualize the porous structure, a three dimensional volume data was built using the set of position and size of bubbles. The simulation box was divided into 10 μm cubes. If a volume of bubble overlaps a cube, the value of the cube become unity otherwise zero. Repeating all bubbles for all cubes in the simulation box, eventually, the volume data is build. Then the volume data is interpolated each other using the smooth3 of MATLAB^® function, which convoluted the volume data with the box kernel. The smoothed volume data was converted to the file format and subjected to a visualization function of ParaView^®. As shown in Fig. 5, the size and number of bubbles in K9 depend on the threshold of contour of volume fraction of matrix. In this study, the threshold of 0.1 was used hereafter.

Fig. 5.

Visualization of volume data of bubbles. The model coal was K9. (Online version in color.)

Moreover, the intensity map of volume fraction of matrix is used to see a cross-section of simulation box. The bubbles in intensity map of Fig. 5 are shown white (red) domains. The numbers of bubbles in the contour and intensity are different because the contour view is the through view of whole simulation box while the intensity view is the view of cross-section.

2.10. Porous Structure of Coke

Two different types of coals were chosen for validation of simulations. K9 coal is a dangerously high coking pressure coal. High MF (Peabody coal) is high-fluidity coal. The characteristics of these coals are shown in Table 2. The prepared coals were charged at a bulk density of 0.90 (g/cm³) (dry basis) in an electrically heated pilot coke oven (430 mm wide, 0.11 m³ internal volume) and carbonized. The detail of preparation of specimen is elsewhere.²⁷⁾

Table 2. Characteristics of coals for observation of porous structure.

Coal	Ash (%) (dry basis)	Volatile matter (%) (dry basis)	Total sulfur (%) (dry basis)	Mean reflectance in oil (%)	Maximum dilatation (%)	Maximum fluidity (ddpm)
K9	9.7	18.5	0.19	1.55	42	51
High MF (Peabody)	6.9	34.4	0.92	0.89	163	20417

3. Results and Discussion

The production of coke is a series of softening and gas generation. To validate our developed model, it is worth to compare coke from two different types of coals: (1) K9 and (3) High MF. The (1) K9 is kept hard (high-viscosity) at elevated temperature and low-gas generation rate coal. The (3) High MF is soft at elevated temperature and produces gas a lot. The comparison of experiment and numerical simulation results of (1) K9 and (3) High MF is necessary to understand the effect of viscosity change as well as the gas generation rate that determines the porous structure of coke.

Figures 6(a-1) and 6(b-1) show cross-sectional pictures of coke of (1) K9 and (3) High MF obtained by the experiments. The cross-section of (1) K9 contains dark 100 μm order domains while that of (3) High MF 500 μm order of dark domains. The dark domains are pores in the coke. (1) K9 is high-viscosity and low gas generation rate so that bubbles’ growth was limited to 100 μm order. (3) High MF is low-viscosity and high-gas generation rate. The bubbles grew sufficiently, coalesced each other and eventually formed heterogeneous porous structure in (3) High MF.

Fig. 6.

Comparison of cross-section of coke pictured in experiments and snapshot of simulation. (a-1): cross-sectional picture of K9 experimentally obtained. (a-2): contour plot of volume fraction of matrix of model coal of K9. (a-3): cross-sectional image of intensity map of model coal of K9. (b-1) cross-sectional picture of High MF experimentally obtained. (b-2): contour plot of volume fraction of matrix of model coal of High MF. (b-3): cross-sectional image of intensity map of model coal of High MF. (Online version in color.)

Snapshots and intensity map of numerical simulations of (1) K9 are shown in Figs. 6(a-2) and 6(a-3), respectively. The spherical bubbles are drawn in the simulation box. The size of bubble of K9 in Figs. 6(a-2) and 6(a-3) is in 100 μm orders. Although the experimental result of K9 is not spherical bubbles, the size of bubbles corresponds to the simulation results of K9. However, the number of pores in experimental result is much higher than that of calculation. The bubbles in the experiment are densely packed than the calculation. It was possible to produce more bubbles by changing the gas generation rate parameter, A in Eq. (16). The bubbles’ growth was delayed by the high-viscosity of (1) K9 while the gas diffusion into the bubbles were increased by the increased gas generation rate. The bubble pressure was increased steeply. Eventually, the calculation of bubble pressure was hardly possible to integrate. Although there was a little space to consider about setting individual parameters of bubble nucleation e.g., f₀ and F for each coals, same values of parameters were used to elucidate effects of viscosity and gas generation rate on porous structure of coke.

The snapshot and intensity map of (3) High MF consist of a number of snowman-like bubbles and chained bubbles. The bubbles were contacted each other during bubble growth period in numerical simulation, they were not coalesced each other and were not formed a unified bubble because the accumulated bubble coalescence time is smaller than that characteristic coalescence time. As the (3) High MF softened up and produced a lot of gas, the number and growth rate of bubbles increased and eventually, final structure of coke became interconnected porous structure.

In our understanding, the calculation and experiment are conformable. There are satellite pores around the large pores. The shape of large pores is not spherical in experiments. In the simulation, the shape of pores is in part of sphere. The mechanical and thermodynamic equilibration induced by interfacial tension and gas diffusion occurs and determines final porous structure in nature. Our developed model, the equilibration is not considered yet because the direct numerical simulation of surfaces should be employed, which raises the calculation cost.

The driving force and impedance of bubble coalescence are the pressure difference and the viscosity, respectively as formulated in Eq. (10). In the direct numerical simulation of bubble coalescence in the manner of hydrodynamics, the time difference for time marching is limited to the ratio of viscosity and pressure. Both high viscosity and large pressure difference at the interface of gas bubble is the bottle neck of calculation. The developed model removed the bottle neck by introducing the characteristic time of bubble coalescence and eased the multi body bubble coalescence simulation. Eventually, the model of bubble coalescence with the visualization technique enabled us to draw the transient structure of bubble coalescence in other words, the chained snowman-like structure. The good coke should be consistent with gas-through-pores and mechanical stiffness. This result of numerical simulation may be used to develop a porous structure for further simulation of transportation of gas through the coke or mechanical testing.

The changes of number density, average diameter and bulk density of bubbles in five different coals are shown in Fig. 7. Figure 8 shows snapshots of coke simulation using five different model coals. The average bubble diameter, number density of bubbles and bulk density are shown in Table 4.

Fig. 7.

Changes of number density, average diameter and bulk density of bubbles in five different coals.

Fig. 8.

Snapshots of coke simulation using five different model coals. The above five images are snapshots of bubbles. The bottom images are the intensity maps of volume fraction of bubbles. The scale bar is 100 μm.

Table 4. Bubble diameter, number density and bulk density of coke simulated using five different model coals.

Coal type	Average bubble diameter [μm]	Number density of bubbles [1/mm3]	Bulk density [g/cm3]
(1) K9	19.5	2237	0.983
(2) Goonyella	98.1	1040	0.561
(3) High MF	103	946	0.543
(4) White Haven	84.4	1310	0.588
(5) ASP	24.4	1972	0.964

The concentration of gas dissolved in coal increased as the function of gas generation rate with temperature. At a certain level of gas concentration at a temperature, a gas bubble is considered to nucleate in the coal as described in Numerical simulation subsection. The curve of number density of bubbles in Fig. 7 started at that temperature. The onset temperatures of bubble nucleation of five different coals are also shown in Fig. 3. The cross point of gas generation curve and the vertical line indicates the gas generation rate at the onset of bubble nucleation. Moreover, the integration of gas generation rate along the temperature divided by the heating rate is the total gas amount generated in the coal. The larger the total gas amount generated is, the lower the onset temperature of bubble nucleation is.

Then the number density of bubbles is accumulated with rising the temperature. Simultaneously, bubbles are expanded by the thermal expansion of gas in the bubbles as well as the diffusion of gas from coal matrix. All of bubble diameters were summed and divided by the total number of bubbles to calculate the average diameter of bubble. The average diameter is increased as the expansion of bubbles. The number density of bubbles is ceased at certain level because the bubble coalescence occurs and bubble nucleation stops. The bubble coalescence and stop of bubble nucleation are induced by the expansion of bubbles, which reduces the distance between the surfaces of bubbles as well as reduces average concentration of gas. The bulk density is decreased by the increase of number density of bubble and the bubble diameters.

The bubble nucleation started in (5) ASP first as the initial gas generation rate is the maximum among five coals. The average bubble size of ASP was limited to 25 μm and the simulation stopped at 423.1°C. The number density of bubbles increased and simultaneously the concentration difference of dissolved gas to bubbles expressed as c–c_R in Eq. (9) was increased by the increases of average gas concentration. At 423.1°C, in calculation, the dissolved gas diffused into bubbles more than net amount of dissolved gas. Therefore, the mass balance of gas was failed. The rapid expansion of bubbles decreased the average concentration while the bubble growth rate was decreased by the decrease of the average concentration. When the mass balance of gas was failed, the deceleration of bubble growth by the decrease of average concentration could not keep up to the consumption of gas to bubble. To perform the simulation further keeping the mass balance of gas, the average concentration should be updated more frequently. The time step of bubble nucleation was needed to reconsider.

Similar results were obtained for (1) K9 and (4) White Haven.

The simulation of (1) K9 was stopped by another reason. The (1) K9 produced the large number of bubbles while the average diameter was relatively small due to high viscosity of (1) K9. With the increasing of average concentration of bubbles as the continuous gas generation and slow bubble growth, the concentration difference between the interface of bubble, c_R, and the average concentration, was increased. The increase of concentration difference increased the change of bubble pressure as described in Eq. (9). Too large change of pressure was difficult to integrate numerically with sufficient accuracy. Eventually, the calculation was failed at the point at 483.2°C.

Although the small size of bubble consumed a little amount of dissolved gas to expand, the total amount of gas became large when the total number of bubbles was fairly large. The drastic decrease of gas occurred by the large number of bubble with a little expansion. Eventually, the total amount of gas in bubbles becomes larger than the total amount of gas produced in coal matrix in the course of simulation.

The minimum viscosity of (4) White Haven was the second highest viscosity among the coal investigated. The gas generation rate of (4) White Haven was also the second highest rate. As shown in Fig. 7, the number density of bubbles increased with the increase of temperature, and then decreased. There was a plateau of number density of bubbles from 433°C to 450°C. On this plateau, the increase of bubbles caused by nucleation was comparable with the decrease of bubbles caused by the bubble coalescence. Because of high-gas generation rate, the bubble nucleation, growth and coalescence are advanced, eventually, the plateau was appeared. After certain decrease of number density caused by the dominant bubble coalescence, the second plateau was appeared. At this plateau, the average bubble diameter also showed a plateau. The bubble nucleation and coalescence was stopped. The bubble growth became dominant, consumed the dissolved gas a lot and finally the mass balance of gas is failed.

As explained above, (1) K9 was the smallest bubble diameter due to high-viscosity and low gas generation rate that caused the lowest total amount of gas. The number density of bubbles was the highest because the generated gas was consumed for bubble nucleation rather than bubble growth of which rate was slow. Comparing (2) Goonyella and (3) High MF, the bubble size of (2) Goonyella was smaller than that of (3) High MF. Both (2) Goonyella and (3) High MF showed many inter-connections between bubbles. The (3) High MF showed inter connections of large and small bubbles. The peak temperature of (3) High MF was 10°C lower than that of (2) Goonyella. As the nucleation occurred faster in (3) High MF, the bubbles could grow longer time than that in (2) Goonyella, which caused to increase the bubble diameter.

There are several problems to perform numerical simulation of K9, White Haven and ASP. These model coals produced huge number of bubbles, which consumed large amount of dissolved gas in the matrix. Eventually, the gas consumption by bubble growth exceeded the amount dissolved in matrix and failed to continue the numerical simulation. Or too extreme increase of bubble pressure occurred to continue the numerical integration of simulation.

In this model, local concentration profile of bubble is overlapped and is disappeared when the two bubbles are coalesced. To account the change of concentration, the average concentration is used to calculate the flux to the bubble. The average concentration makes the uneven concentration in matrix flat. This approximation is too simple to conduct the simulations for huge bubble nucleation cases.

Moreover, the rule-based model of bubble coalescence still has a problem where new position of two coalesced bubbles. In the model, it cannot avoid that the new position is possibly located inside of bubble. Two bubble shares a part of volume of their bubbles.

Furthermore, this model assumed that the generated gas was methane. However, there are several kinds of gas e.g., CO, CO₂ and H₂ which should be involved by kinetic model of coal pyrolysis reactions.

4. Conclusion

A simplified numerical simulation model of bubble nucleation, growth and coalescence was developed for understanding the coke production process from coal. The bubble nucleation caused by gas generation in coke, bubble growth induced by gas-diffusion, thermal expansion of gas and pressure difference across the bubble, bubble coalescence between two bubbles were mathematically expressed and solved numerically in the course of elevating temperature in coal. To model the bubble coalescence of two bubbles, a rule-based model where the coalescence depended on the viscosity and pressure difference of two bubbles was developed. The rule could be embedded in the simulation to decide which pair of bubbles was likely coalesced.

In spite of simplified model, the developed model resulted in snapshots of porous structure of coke using five different model coals e.g., K9, Goonyella, High MF, White Haven and ASP. The cross-sectional pictures of Goonyella and High MF obtained by experiments were corresponded to the snapshots of simulation qualitatively in the view points of inter-connectivity of bubbles and their size.

Nomenclature

A: constant of gas generation rate

D: diffusion coefficient

D_c: conversion factor of viscosity to diffusion coefficient.

F: correction factor of free energy barrier of bubble nucleation

J: bubble nucleation rate

N_A: the Avogadro number

$P_{\mathrm{C}}^{*}$ : pressure outside of nucleus gas bubble

P_D: bubble pressure

P_C: pressure outside of gas bubble.

$P_{\mathrm{D}}^{*}$ : mechanically equilibrium pressure of gas bubble

P_D0: bubble pressure of nucleus

P_D1: bubble pressure of bubble 1

P_D2: bubble pressure of bubble 2

R: radius of bubble

R₀: bubble radius of nucleus

R_i: bubble radius of bubble i

T: absolute temperature

T_C: peak temperature

c: concentration of gas dissolved in matrix

c_∞: gas concentration in infinity

c_R: gas concentration in the matrix at the surface of bubble

c₀: initial concentration of gas in matrix

c: average gas concentration dissolved in matrix

f₀: correction factor of probability of bubble nucleation

k_B: the Boltzmann constant

k_H: solubility parameter of gas

m: mass of gas molecule

n_g: number mole of gas molecules in bubble

$\Re$ : gas constant

t: time

v₀: volume of matrix

w: temperature width of half of peak.

x_i: position vector of bubble i

y(t): gas generation rate

γ: surface tension of polymer

π: circular constant

η: viscosity

δ: width of concentration profile

References

• 1)   D.  Hays,  J. W.  Patrick and  A.  Walker: Fuel, 55 (1976).
• 2)   H.  Aoki,  K.  Goto,  M.  Ueki,  T.  Miura,  K.  Fukuda and  S.  Itagaki: Tetsu-to-Hagené, 82 (1996), 404.
• 3)   Y.  Takizawa,  T.  Kimura,  K.  Higuchi,  Y.  Morozumi,  H.  Aoki,  T.  Miura,  K.  Uebo and  K.  Fukuda: Testu-to-Hagané, 92 (2006), 198.
• 4)   M. S.  Oh,  W. A.  Peters and  J. B.  Howard: AIChE J., 35 (1989).
• 5)   K.  Taki,  D.  Kitano and  M.  Ohshima: Ind. Eng. Chem. Res., 50 (2011), 3247.
• 6)   K.  Taki and  S.  Okumura: Macromolecules, 43 (2010), 9899.
• 7)   K.  Taki: Chem. Eng. Sci., 63 (2008), 3643.
• 8)   K.  Taki,  K.  Tabata,  S.  Kihara and  M.  Ohshima: Polym. Eng. Sci., 46 (2006), 680.
• 9)   K.  Taki,  K.  Nitta,  S.-I.  Kihara and  M.  Ohshima: J. Appl. Polym. Sci., 97 (2005), 1899.
• 10)   K.  Taki,  T.  Yanagimoto,  E.  Funami,  M.  Okamoto and  M.  Ohshima: Polym. Eng. Sci., 44 (2004), 1004.
• 11)   K.  Taki,  T.  Nakayama,  T.  Yatsuzuka and  M.  Ohshima: J. Cell. Plast., 39 (2003), 155.
• 12)   A.  Toramaru: J. Geophys. Res. [Atmospheres], 100 (1995), 1913.
• 13)   M.  Blander and  J.  Katz: AIChE J., 21 (1975), 833.
• 14)   D. E.  Rosner and  M.  Epstein: Chem. Eng. Sci., 27 (1972), 69.
• 15)   P.  Payvar: Int. J. Heat Mass Tran., 30 (1987), 699.
• 16)   C.  Han and  H.  Yoo: Polym. Eng. Sci., 21 (1981), 518.
• 17)   M. A.  Shafi,  K.  Joshi and  R. W.  Flumerfelt: Chem. Eng. Sci., 52 (1997), 635.
• 18)   R. B.  Bird,  R. C.  Armstrong and  O.  Hassager: Dynamics of Polymeric Liquids, Vol. 1, Fluid Mechanics, John Wiley & Sons, Inc., New Jersey, (1987).
• 19)   K.  Miura and  T.  Maki: J. Chem. Eng. Jpn., 31 (1998), 228.
• 20)   J.  Hayashi,  M.  Morita,  H.  Kumagai and  T.  Chiba: Fuel, 82 (2003), 487.
• 21)   A. K.  Doolittle: J. Appl. Phys., 22 (1951), 1471.
• 22)   M. L.  Williams,  R. F.  Landel and  J. D.  Ferry: J. Am. Chem. Soc., 77 (1955), 3701.
• 23)   J. S.  Vrentas and  J. L.  Duda: J. Polym. Sci. B, 15 (1977), 403.
• 24)   J. S.  Vrentas and  J. L.  Duda: J. Polym. Sci. B, 15 (1977), 417.
• 25)   J. S.  Vrentas and  J. L.  Duda: J. Polym. Sci. B, 15 (1977), 441.
• 26)   Y.  Sato,  K.  Fujiwara,  T.  Takikawa,  Sumarno,  S. Takishima and H. Masuoka: Fluid Phase Equilibr., 162 (1999), 261.
• 27)   T.  Arima: Testu-to-Hagané, 87 (2001), 274.