2.1. Pore-scale domain with realistic coke deposition pattern

Coke combustion was investigated in a two-dimensional (2-D) pore-scale computational domain constructed from micro-tomographic images by micro-CT (nanoVoxel-3000) for the coked porous media. The micro-CT system built digital images of 1700 × 1700 × 2100 voxels with a voxel resolution of 500 nm. A 2-D slice consisting of 1200 × 2000 voxels (dimensions 0.6 × 1 mm²) or 2.4 million grid cells, shown in figure 2(a), was cropped from the reconstructed data and imported into the ilastik toolkit (Berg et al. Reference Berg2019) for image segmentation by machine learning algorithms. The image segmentation workflow began from label assignments of the coke, rock and pore pixels through human annotations drawn on the raw greyscale image. Random forest classifiers were trained to learn the labelled features and empower the prediction capacity to partition the image into the coke, rock and pore phases. The labelling–training–prediction workflow was looped in an interactive mode to correct the classification flaw until the desired segmentation was achieved, as illustrated in figure 2(b). After segmentation, the volume fraction of the pore pixels was 23.8 %, while the coke volume fraction was 8.19 %. For the present coke combustion simulation, two additional open flow buffer regions 20 pixels wide were added before and after the entrance of the coked porous medium so that the fully developed boundary condition could be enforced at the inlet and outlet boundaries of the computational domain.

Figure 2. (a) Two-dimensional raw image slice scanned by micro-CT showing rock grains (lightest grey), coke (grey) and pore (darkest grey) and (b) segmented image slice for computational domain with coloured phase of interest for rock grains (blue), coke (red) and pore (black).

As observed in figure 2(b), the base matrix was packed with an idealized sphere-like glass bead pack, regarded as the sandstone proxy. Coke tended to deposit inside the pore throat, which agrees with previous observations from scanning electron microscopy (Xu et al. Reference Xu, Long, Jiang, Ma, Zan, Ma and Shi2018a). The coke deposition pattern blocked all of the connectable pore network if the coke was considered an impermeable obstacle. This implied that the previous numerical method (Xu et al. Reference Xu, Long, Jiang, Ma, Zan, Ma and Shi2018a), which neglected the effect of the sub-resolution continuum on the fluid flow, is not suitable for modelling the reactive flow in such a realistic coked porous medium. Additionally, the non-uniform spatial distribution of the coke deposition resulted in coke accumulations that were rich in some local regions but poor in other regions, increasing the heterogeneity of the coked porous medium. The naturally heterogeneous coke deposition in pore throats showed a significant difference compared with the previous idealized structures (Xu et al. Reference Xu, Long, Jiang, Zan, Huang, Chen and Shi2018b; Lei et al. Reference Lei, Wang and Luo2021), which were composed of staggered or random arrangements of circle-like grains with all coke deposited on the grain surface.

2.2. Governing equation and assumptions

The micro-CT system with micrometre-scale resolution and millimetre-scale field of view cannot resolve the pore structure in coke regions because the characteristic length scales of the coke nanopores and the rock macropores differ by orders of magnitudes. In addition, the direct pore-scale simulation for all the scales is unacceptable owing to the enormous computational cost for billion-cell representations despite the 2-D simulation problems. Therefore, the multiscale modelling framework (i.e. the micro-continuum approach for pore-scale simulation) was introduced to establish the governing equations for the multiscale coke combustion dynamics. With the micro-continuum approach, structured Cartesian-type grid cells directly correspond to the voxel elements of the segmented image in figure 2(b) without requiring a complex mesh technique. The volume fraction of void space in each grid cell, ${\varepsilon _f}$, referred to as the local porosity, was imposed to represent the complex multiscale pore structure and filter the explicit structural information below the micro-CT resolution limit. In the current filtering strategy, the pore regions labelled in the segmented image were denoted by ${\varepsilon _f} = 1$ to represent the clear fluid zone, while the rock regions were denoted by ${\varepsilon _f} = 0$ to represent impermeable solid zones without small pore space. For the coke regions, an intermediate value of $0 < {\varepsilon _f} < 1$ can be given to symbolize the sub-resolution porous medium with an aggregate of nanopores and coke solids, as shown in figure 1. The intermediate value is equal to the intrinsic coke porosity, which depends on coke formation conditions, such as crude oil properties and pyrolysis temperature. In the present study, the sub-resolution coke porosity was specified as 0.7 to account for the gas transport and chemical reactions in such porous coke structures. The grid cells were assumed to not contain both the rock and the coke; therefore, the volume fraction of the rock phase, ${\varepsilon _{rock}}$, or the volume fraction of the coke phase, ${\varepsilon _{coke}}$, can be derived according to the known values of ${\varepsilon _f}$ and the constraint relation as

(2.1)\begin{equation}{\varepsilon _f} + {\varepsilon _{coke}} + {\varepsilon _{rock}} = 1.\end{equation}

In the micro-continuum approach, partial differential equations (PDEs) are averaged over control volumes, corresponding to the grid cells, to govern the conservation law of energy, mass and species concentrations, which are modelled with averaged field variables. According to the volume averaging theory, these PDEs hold for all computational cells regardless of their contents. Before deriving the governing equations, the following assumptions were made to simplify the modelling of the complex hybrid-scale physics:

1. (1) the non-slip boundary condition is held at the impermeable rock interface owing to the approximated Knudsen number of 2.28 × 10⁻³ in the resolved macropores surrounded by the rock grains;

2. (2) the heat capacities of the gas and solid components are constant and only depend on the initial temperature and pressure (Lei et al. Reference Lei, Wang and Luo2021);

3. (3) the densities of the solid components are constant;

4. (4) the gas is an ideal gas with the equation of state as (2.2)

   (2.2)\begin{equation}{\rho _f} = {\bar{p}_f}/R\bar{T},\end{equation}

where ${\bar{p}_f}$ is the intrinsically phase averaged pressure, ${\bar{p}_f} = (1/{V_f})\int_{{V_f}} {{p_f}} \,\textrm{d}V$;

1. (5) the radiative heat transfer is not considered;

2. (6) the thermal equilibrium condition is valid within each coke cell at 500 nm resolution when the burning temperature is less than 1200 K, which maintains

   (2.3)\begin{equation}\varUpsilon = \frac{{\left|{\frac{{\dot{r}{h_r}}}{{\rho {c_p}}}} \right|}}{{\left|{\alpha \frac{{{\partial^2}T}}{{\partial {x^2}}}} \right|}} = O\left( {\frac{{\dot{r}{h_r}{L^2}}}{{\rho {c_p}\alpha T}}} \right) \ll 1,\end{equation}

where $\varUpsilon$ is the ratio of the reaction term and the thermal conductive term in the energy equation, h_r is the coke combustion heat, $\dot{r}$ is the coke combustion rate, and L is the cell size, 500 nm;

1. (7) the coke combustion is modelled as (2.4), which indicates that the coke reacts with the oxygen in the injected hot air to produce CO₂. The reaction rate is estimated according to a first-order Arrhenius equation as (2.5)

   (2.4)\begin{gather}\textrm{C} + {\textrm{O}_\textrm{2}} \to \textrm{C}{\textrm{O}_\textrm{2}} + {h_r},\end{gather}
   (2.5)\begin{gather}\dot{r} = {a_v}Aexp \left( { - \frac{E}{{R\bar{T}}}} \right)\frac{{{\rho _f}{{\bar{\omega }}_{f\textrm{,}{\textrm{O}_\textrm{2}}}}}}{{{M_{{\textrm{O}_\textrm{2}}}}}},\end{gather}

where ${a_v}$ is the reactive surface area per control volume, which is only non-zero at the grid cells containing the coke. The detailed sub-grid model of the reactive surface area will be discussed in § 2.3. In the coke combustion kinetics formulation, A and E are the pre-exponential factor and activation energy, respectively, while R is the ideal gas constant, ${\rho _f}$ is the gas density, and ${M_{{\textrm{O}_\textrm{2}}}}$ is the molecular weight of the oxygen. To match the following volume averaging governing equations, $\bar{T}$ represents the temperature averaged over the control volume, $\bar{T} = (1/V)\int_V T \,\textrm{d}V$, and ${\bar{\omega }_{f\textrm{,}{\textrm{O}_\textrm{2}}}}$ is the O₂ mass fraction intrinsically averaged over the gas phase, ${\bar{\omega }_{f\textrm{,}{\textrm{O}_\textrm{2}}}} = (1/{V_f})\int_{{V_f}} {{\omega _{f\textrm{,}{\textrm{O}_\textrm{2}}}}} \,\textrm{d}V$. In this simulation, the coke combustion kinetics of A, E and h_r are set as $9.717 \times {10^6}\;\textrm{m}\;{\textrm{s}^{ - 1}}$, 131.09 kJ mol⁻¹ and 388.5 kJ mol⁻¹, respectively (Ren, Freitag & Mahinpey Reference Ren, Freitag and Mahinpey2007), which are the same values as in previous studies (Xu et al. Reference Xu, Long, Jiang, Zan, Huang, Chen and Shi2018b; Lei et al. Reference Lei, Wang and Luo2021) for better comparisons of different numerical models.

According to the above analysis and assumptions, a set of the volume averaging governing equations for multiscale coke combustion are shown in (2.6)–(2.10), which are valid throughout the domain regardless of the content of the grid cell:

(2.6)\begin{gather}\frac{1}{{{\varepsilon _f}}}\left( {\frac{{\partial {\rho_f}\bar{\boldsymbol{U}}}}{{\partial t}}+ \boldsymbol{\nabla }\boldsymbol{\cdot }\left( {\frac{{{\rho_f}}}{{{\varepsilon_f}}}\boldsymbol{\bar{U}\bar{U}}} \right)} \right) ={-} \boldsymbol{\nabla }{\bar{p}_f} + {\rho _f}g + \frac{{{\mu _f}}}{{{\varepsilon _f}}}{\nabla ^2}\bar{\boldsymbol{U}} - {\mu _f}{k^{ - 1}}\bar{\boldsymbol{U}},\end{gather}

(2.7)\begin{gather}\frac{{\partial {\rho _f}{\varepsilon _f}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\rho _f}\bar{\boldsymbol{U}}) ={-} {\dot{m}_{coke}},\end{gather}

(2.8)\begin{gather}\frac{{\partial {\rho _{coke}}{\varepsilon _{coke}}}}{{\partial t}} = {\dot{m}_{coke}},\end{gather}

(2.9)\begin{gather}\frac{{\partial {\varepsilon _f}{\rho _f}{{\bar{\omega }}_{f,i}}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\rho _f}\boldsymbol{U}{\bar{\omega }_{f,i}}) = \boldsymbol{\nabla }\boldsymbol{\cdot }({\varepsilon _f}{\rho _f}{\hat{D}_{i,eff}}\boldsymbol{\nabla }{\bar{\omega }_{f,i}}) + {\dot{m}_i},\quad i = {\textrm{O}_\textrm{2}}\textrm{,C}{\textrm{O}_\textrm{2}},\end{gather}

(2.10) \begin{gather}\dfrac{{\partial {\varepsilon _s}{\rho _s}{h_s}}}{{\partial t}} + \dfrac{{\partial {\varepsilon _f}{\rho _f}{h_f}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\rho _f}\bar{\boldsymbol{U}}{h_f}) + \dfrac{{\partial {\rho _f}K}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\rho _f}\bar{\boldsymbol{U}}K)\nonumber\\ \qquad= \dfrac{{\partial {{\bar{p}}_f}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\lambda _{eff}}\boldsymbol{\nabla }T) + {\rho _f}\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{g} + {{\dot{q}}_r}.\end{gather}

The momentum conservation equation, known as the DBS equation (Brinkman Reference Brinkman1949a), is expressed in (2.6), where $\bar{\boldsymbol{U}}$ is the superficial velocity, namely, the Darcy velocity, t is the time, g is the gravitational acceleration and ${\mu _f}$ is the gas viscosity. The k ⁻¹ in the last term of the (2.6) is the reciprocal of the absolute permeability, which varies with the local porosity via the commonly used Kozeny–Carman equation:

(2.11)\begin{equation}{k^{ - 1}} = k_0^{ - 1}\frac{{{{(1 - {\varepsilon _f})}^2}}}{{\varepsilon _f^3}}.\end{equation}

According to the Kozeny–Carman relation, the porous resistance term, ${\mu _f}{k^{ - 1}}\bar{\boldsymbol{U}}$, vanishes in the resolved macropore region (${\varepsilon _f} = 1,\;{k^{ - 1}} = 0\;{\textrm{m}^{ - \textrm{2}}}$) so that the DBS equation degenerates to the standard Navier–Stokes equation to model the weakly compressible gas flow in the pore space. In the rock regions, a very small decimal number ${\varepsilon _f} = 0.01$ was used to replace ${\varepsilon _f} = 0$ to escape floating-point exceptions, leading to the permeability k of $1.02 \times {10^{ - 3}}\;\textrm{mD}$ with the specified parameter ${k_0}$ of 10⁻¹² in the Kozeny–Carman relation. The rock permeability of $1.02 \times {10^{ - 3}}\;\textrm{mD}$ is consistent with the properties of the impermeable rock in the literature (Lis-Śledziona Reference Lis-Śledziona2019), which reported the lowest permeability of 10⁻³ mD in their interpreted wells. More sensitivity analysis of ${\varepsilon _f}$ on the coke combustion dynamics and the numerical stability can be found in the supplementary material. The velocity in the rock zones reduces to almost zero so that the DBS equation can recover both the impermeable properties of the internal rock regions and the non-slip boundary conditions at the fluid–rock interface (Soulaine & Tchelepi Reference Soulaine and Tchelepi2016). For the porous coke regions, the DBS equation is reduced towards the Darcy law because the porous resistance term gradually becomes dominant over the dissipative viscous term with $0 < {\varepsilon _f} < 1$ (Tam Reference Tam1969; Soulaine & Tchelepi Reference Soulaine and Tchelepi2016). When the coke is about to burn out, only a few fragments remain in the sub-resolution continuum, which leads to a porosity of approximately one. Under such conditions, the investigation (Lévy Reference Lévy1983; Auriault Reference Auriault2008) showed that the flow law obeys the Darcy–Brinkmann law with an approximation of $O(\varepsilon )$. We also recognized that the Kozeny–Carman equation might not be the best choice to relate the permeability and the porosity for coke nanoporous structures. With respect to the problem under consideration, however, the gas transport by Knudsen diffusion was more significant than that by viscous flow, which was located in the transition regime with a Knudsen number of $O(1)$. Therefore, the Kozeny–Carman equation was employed to assign an intrinsic permeability of $O(1)\sim O(10)\;\textrm{mD}$ to the coke regions for weak convective fluxes. Therefore, the DBS equation allows for the Darcy/Darcy–Brinkmann and Navier–Stokes flow to be simulated in the present multiscale coked porous medium (Brinkman Reference Brinkman1949b; Neale & Nader Reference Neale and Nader1974; Whitaker Reference Whitaker1986; Golfier et al. Reference Golfier, Zarcone, Bazin, Lenormand, Lasseux and Quintard2002, Reference Golfier, Bazin, Lenormand and Quintard2004; Rein Reference Rein2009; Luo et al. Reference Luo, Quintard, Debenest and Laouafa2012, Reference Luo, Laouafa, Guo and Quintard2014, Reference Luo, Laouafa, Debenest and Quintard2015; Scheibe et al. Reference Scheibe, Perkins, Richmond, McKinley, Romero-Gomez, Oostrom, Wietsma, Serkowski and Zachara2015; Soulaine & Tchelepi Reference Soulaine and Tchelepi2016).

The mass conversion equations for the gas and coke phases are written as (2.7) and (2.8), respectively, where ${\rho _{coke}}$ is the coke density and ${\dot{m}_{coke}}$ is the coke burning rate, computed by ${\dot{m}_{coke}} ={-} {M_{coke}}\dot{r}$. Here, ${\dot{m}_{coke}}$ is a negative value because the solid coke gradually burns with O₂ to release CO₂ into the gas phase as the volume fraction of the coke ${\varepsilon _{coke}}$ decreases. Therefore, the coke morphology continually evolves with coke combustion, while the mass of the gas phase progressively increases.

The evolution of the species mass fraction, including O₂ and CO₂, is described by the volume averaging convection–diffusion equation in (2.9) with some additional insignificant terms being discarded (Golfier et al. Reference Golfier, Zarcone, Bazin, Lenormand, Lasseux and Quintard2002), where ${\bar{\omega }_{f,i}}$ is the mass fraction of species components in gas and ${\hat{D}_{i,eff}}$ is the effective diffusivity of species i. Here, ${\hat{D}_{i,eff}}$ is given by

(2.12) \begin{equation}{\hat{D}_{i,eff}} = \left\{ {\begin{array}{*{20}{@{\quad}l}} {{{\hat{D}}_{i,b}},}&{\textrm{in}\;\textrm{the}\;\textrm{resolved}\;\textrm{macropore}\;\textrm{regions}}\\ {\dfrac{{{\varepsilon_f}}}{\tau } \times \dfrac{1}{{\left( {\dfrac{1}{{{{\hat{D}}_{i,b}}}} + \dfrac{1}{{{{\hat{D}}_{i,Knud}}}}} \right)}}\textrm{, }}&{\textrm{in}\;\textrm{the}\;\textrm{unresolved}\;\textrm{nanopore}\;\textrm{regions}} \end{array}} \right.\end{equation}

where ${\hat{D}_{i,b}}$ is the effective diffusion coefficient for component i of the mixture in the resolved macropore regions, ${\hat{D}_{i,Knud}}$ is the Knudsen diffusion coefficient and $\tau$ is the tortuosity of the coke nanoporous structure, taken as 1 owing to the cylindrical pores in coke. The effective diffusion coefficient ${\hat{D}_{i,b}}$ can be computed by the Wilke relation (Fairbanks & Wilke Reference Fairbanks and Wilke1950) as (2.13). The average pore radius in the sub-resolution nanopores is of the order of magnitude of 10 nm (Fei et al. Reference Fei, Hu, Xiang, Sun, Fu and Chen2011; Wang et al. Reference Wang, Liu, Li, Yang, Wang and Song2018), which leads to a Knudsen number of $O(1)$ and a transition regime. Therefore, Knudsen diffusion should be considered for the coke regions and can be determined by kinetic theory as (2.14). The well-known Bosanquet relation (Krishna & Wesselingh Reference Krishna and Wesselingh1997), as the second term of (2.12), was introduced to combine the effect of molecular diffusion and Knudsen diffusion, which is consistent with the dusty-gas model with acceptable deviations (Veldsink et al. Reference Veldsink, Van Damme, Versteeg and Van Swaaij1995). The additional factor ${\varepsilon _f}/\tau$ before the Bosanquet formula accounts for diffusion resistances owing to the porous structures (Wakao & Smith Reference Wakao and Smith1962; Soulaine & Tchelepi Reference Soulaine and Tchelepi2016; Soulaine et al. Reference Soulaine, Roman, Kovscek and Tchelepi2017). The mechanical dispersion in the sub-resolution continuum was not coupled in the diffusion term owing to its insignificance compared with the molecular and Knudsen diffusion under the low Péclet number during coke combustion (Bijeljic & Blunt Reference Bijeljic and Blunt2007),

(2.13)\begin{gather}{\hat{D}_{i,b}} = \frac{{1 - {{\bar{\chi }}_{f,i}}}}{{\sum\limits_{j{\ne}i} {({{\bar{\chi }}_{f,j}}/{{\hat{D}}_{ij}})} }},\quad i = {\textrm{O}_\textrm{2}}\textrm{,C}{\textrm{O}_\textrm{2}}\textrm{,}{\textrm{N}_\textrm{2}},\end{gather}

(2.14)\begin{gather}{\hat{D}_{i,Knud}} = \frac{{2{d_p}}}{3}\sqrt {\frac{{2R\bar{T}}}{{\mathrm{\pi }M}}} ,\end{gather}

where ${\bar{\chi }_{f,i}}$ is the molar fraction of component i, ${\hat{D}_{ij}}$ is the binary diffusion coefficient of component i in component j (Cussler & Cussler Reference Cussler and Cussler2009) and ${d_p}$ is the average pore diameter, specified as 20 nm (Fei et al. Reference Fei, Hu, Xiang, Sun, Fu and Chen2011; Wang et al. Reference Wang, Liu, Li, Yang, Wang and Song2018). The source term in (2.9) represents the O₂ consumption rate or CO₂ production rate owing to coke combustion, which is given by ${\dot{m}_{{\textrm{O}_\textrm{2}}}} ={-} {M_{{\textrm{O}_\textrm{2}}}}\dot{r}$ or ${\dot{m}_{\textrm{C}{\textrm{O}_\textrm{2}}}} = {M_{\textrm{C}{\textrm{O}_\textrm{2}}}}\dot{r}$. Similar to the DBS equation, (2.9) tends asymptotically to the ordinary advection–diffusion equation in the resolved macropores but also considers the effect of gas rarefaction in the unresolved nanopores without separate solvers and domains.

The thermal equilibrium condition is assumed to hold within each control volume, which consists of solid coke and its internal micropore voids. This allows the energy conservation equations for the solid and fluid phase to be combined into a single equation, as shown in (2.10), where ${h_s}$ and ${h_f}$ are the specific enthalpy for the solid phase and fluid phase, respectively, K is the specific kinetic energy, $K = |\bar{\boldsymbol{U}}{|^2}/2$, and ${\lambda _{eff}}$ is the effective thermal conductivity, given by ${\lambda _{eff}} = {\varepsilon _f}{\lambda _f} + {\varepsilon _{coke}}{\lambda _{coke}} + {\varepsilon _{rock}}{\lambda _{rock}}$ (Soulaine & Tchelepi Reference Soulaine and Tchelepi2016). The last term of (2.10), ${\dot{q}_r}$, corresponds to the heat source term released by exothermic coke combustion, which can be calculated by ${\dot{q}_r} = {h_r}\dot{r}$. Equation (2.10) cannot be directly solved because it contains two unknown field variables, ${h_s}$ and ${h_f}$. Along with the constant heat capacity assumption, ${h_s}$ can be obtained by ${h_s} = {C_{p,s}}T = {C_{p,s}}({h_f}/{C_{p,f}})$; therefore, the energy balance equation can be recast into (2.15), which retains ${h_f}$ as the primary variable. In the simulation, the thermophysical properties of the gas species, such as the specific heat capacity, viscosity and Prandtl number, came from REFPROP Version 8.0 (Lemmon, Huber & McLinden Reference Lemmon, Huber and McLinden2002), while the thermophysical properties of the rock phase were the same as in Xu et al. (Reference Xu, Long, Jiang, Zan, Huang, Chen and Shi2018b),

(2.15) \begin{gather} \dfrac{{\partial {\varepsilon _f}{\rho _f}{h_f}}}{{\partial t}} + \dfrac{{\partial {\varepsilon _{coke}}({\rho _{coke}}{C_{p,coke}}/{C_{p,f}}){h_f}}}{{\partial t}} + \dfrac{{\partial {\varepsilon _{rock}}({\rho _{rock}}{C_{p,rock}}/{C_{p,f}}){h_f}}}{{\partial t}}\nonumber\\ \quad + \,\boldsymbol{\nabla }\boldsymbol{\cdot }({\rho _f}\bar{\boldsymbol{U}}{h_f}) + \dfrac{{\partial {\rho _f}K}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\rho _f}\bar{\boldsymbol{U}}K)\nonumber\\ \qquad\quad = \dfrac{{\partial {{\bar{p}}_f}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left( {\textrm{ }\dfrac{{{\lambda_{eff}}}}{{{C_{p,f}}}}\boldsymbol{\nabla }{h_f}} \right) + {\rho _f}\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{g} + \dot{q}.\end{gather}

For boundary conditions, as illustrated in figure 2(b), the inlet is imposed by the Poiseuille velocity profile and the Dirichlet boundary condition for the temperature and mass fraction of gas species, which contain a mixture of 78 % N₂ and 22 % O₂. At the outlet of the computational domain, fully developed thermal and reactive flows are enforced. At the top and bottom boundaries, the wall-type boundary is specified with adiabatic conditions.

The developed pore-scale micro-continuum model involves more physics than the previous LB models, including modelling of the sub-resolution reactive transport, two-way coupling between the chemical reaction and fluid flow, improved energy conservation during the structural evolution, and fewer assumptions in the thermal and transport properties. During the structural evolution of the sub-resolution continuum, the effective medium coefficients were also temporally updated, including the permeability, effective diffusivity and effective thermal conductivity. More discussions on these improvements can be found in the supplementary material. The current numerical model enforced the no-slip velocity condition on the interfaces between the resolved macropores and the rock grains for the high-permeability sandstone-like rock under consideration. However, when the numerical model is applied to simulate coke combustion through low-permeability carbonate rock or under low-pressure conditions with increasing Knudsen number, the slip or even transition flow regime should be accounted for in the resolved mesopores and possibly nanopores so that the slip velocity condition can be recovered at the rock walls. The unified pore-scale micro-continuum approach that applies to all flow regimes can be found in these studies (Guo, Ma & Tchelepi Reference Guo, Ma and Tchelepi2018; Soulaine et al. Reference Soulaine, Creux and Tchelepi2019).

2.3. Reactive surface area model

The variables ${\varepsilon _f}$, k and ${\hat{D}_{i,eff}}$ in the volume averaging governing models above are used to model the flow and diffusion behaviours within the coke cells with sub-resolution nanopores. The sub-grid models of the permeability k and the diffusion coefficient ${\hat{D}_{i,eff}}$ are introduced in § 2.2 to quantify these transport properties with respect to the local porosity, ${\varepsilon _f}$. The last component to complete the mathematical models is the sub-grid model of the reactive surface area within the coke. The integration of the reactive surface area model with micro-continuum models can avoid the explicit treatment of the reactive fluid–solid interface with complex interfacial conditions. In the present study, the diffuse interface model (DIM) (Luo et al. Reference Luo, Quintard, Debenest and Laouafa2012) and the random pore model (RPM) (Bhatia & Perlmutter Reference Bhatia and Perlmutter1980) are presented to investigate the effect of the reactive surface area models on the coke combustion dynamics. The diffuse interface model is commonly adopted in the mineral dissolution problem to reproduce a sharp diffuse interface where the dissolution phenomenon occurs (Luo et al. Reference Luo, Quintard, Debenest and Laouafa2012; Soulaine et al. Reference Soulaine, Roman, Kovscek and Tchelepi2017). In this work, the diffuse interface model is expressed as

(2.16)\begin{equation}{a_v} = 2||\boldsymbol{\nabla }{\varepsilon _{coke}}||(1 - \varepsilon _f^2),\end{equation}

where the term $2||\boldsymbol{\nabla }{\varepsilon _{coke}}||$ is the gradient of the coke volume fraction with a correction factor of 2, which characterizes the average interfacial surface of a control volume based on the volume averaging theorem (Whitaker Reference Whitaker2013). As a result, reactive surface areas are non-zero only at the gas–coke interface. The term $(1 - \varepsilon _f^2)$ controls the diffuse thickness, which follows the recommended formulation in Luo's work (Luo et al. Reference Luo, Quintard, Debenest and Laouafa2012). Therefore, the diffuse interface model can serve as an immersed boundary condition to recover the reactive flux boundary condition at the edge of the coke regions.

The random pore model is written as (Bhatia & Perlmutter Reference Bhatia and Perlmutter1980)

(2.17)\begin{equation}{a_v} = (1 - x){a_v}_0\sqrt {1 - \psi \ln (1 - x)} ,\end{equation}

where x denotes the coke conversion rate, ${a_v}_0$ is the initial pore surface area of coke and $\psi$ is the called structural parameter, set as 10 in this work. The random pore model assumes that the porous coke structure comprises cylindrical pores with arbitrary pore size distribution and orientations. In contrast to the diffuse interface model, the heterogeneous gas–solid reaction is not limited to the edge of the coke cell but occurs on the internal surfaces of cylindrical nanopores. Therefore, the effective surface areas derived from the random pore model are non-zero within the coke cells. From (2.17), the effective surface area changes as a function of the coke conversion rate, x. The random pore model is formulated to model the initial growth of the reaction surface with increasing conversion as experimentally observed (Bhatia & Perlmutter Reference Bhatia and Perlmutter1980). The maximum reaction surface occurs at a conversion rate of 0.393 (Bhatia & Perlmutter Reference Bhatia and Perlmutter1980), while the enlargement ratio relative to the initial pore surface area can be adjusted by the structural parameter $\psi$.

3.1. Equation discretization

The solver was developed based on the open-source CFD software OpenFOAM (Jasak Reference Jasak1996; Weller et al. Reference Weller, Tabor, Jasak and Fureby1998; Jasak, Jemcov & Tukovic Reference Jasak, Jemcov and Tukovic2007) to solve the governing equations and the sub-grid models formed by (2.6)–(2.17) with the finite-volume method (FVM). In the FVM, the governing equations were first discretized by integrating over each control volume at the time step to yield a set of algebraic equations. The first-order Euler time scheme was used to discretize time derivative $\partial /\partial t$ terms, while the spatial terms were discretized in second-order numerical schemes. The gradient $\boldsymbol{\nabla }$ term was assigned the Gauss linear scheme with a linear scheme implemented for the value interpolation from cell centres to face centres. The Gauss linearUpwind scheme was employed for the divergence $\boldsymbol{\nabla }\boldsymbol{\cdot }(({\rho _f}/{\varepsilon _f})\boldsymbol{\bar{U}\bar{U}})$ term, while other divergence terms for advection of scalar fields, such as $\boldsymbol{\nabla }\boldsymbol{\cdot }({\rho _f}\bar{\boldsymbol{U}}{h_f})$, were discretized with the bounded Gauss limitedLinear scheme to ensure the boundedness. The Laplacian term for thermal diffusion $\boldsymbol{\nabla }\boldsymbol{\cdot }(({\lambda _{eff}}/{C_{p,f}})\boldsymbol{\nabla }{h_f})$ was calculated using the Gauss harmonic orthogonal scheme to guarantee diffusion flux continuity at the fluid–solid interface during conjugate heat transfer, while the Gauss linear orthogonal scheme was used for other Laplacian terms, such as $\boldsymbol{\nabla }\boldsymbol{\cdot }({\varepsilon _f}{\rho _f}{\hat{D}_{i,eff}}\boldsymbol{\nabla }{\bar{\omega }_{f,i}})$.

3.2. Solution algorithms

Sequential coupling strategies were introduced to solve the discretized governing equations for the nonlinear problem, as illustrated in figure 3. In the numerical workflow, the pressure–velocity coupling in the transient non-isothermal reactive flow was handled by the iterative PIMPLE algorithm to improve numerical stabilities, which combined the pressure implicit with splitting of operator (PISO) algorithm (Issa Reference Issa1986) and semi-implicit method for pressure-linked equations (SIMPLE) algorithm (Patankar & Spalding Reference Patankar and Spalding1983). The equation discretization for the PIMPLE algorithm can be found in Appendix A.

Figure 3. Flow chart of the numerical models for the coupled thermal and reactive transport.

The main iterative procedure to advance the time step from n to n + 1 is described as follows:

1. (1) the coke reaction rate ${\dot{r}^\ast }$ is calculated by (2.5) with the temperature, mass fraction of O₂ and specific surface area from the previous iteration step or time step;

2. (2) the volume fraction of coke, $\varepsilon _{coke}^\ast $, is explicitly solved by (2.8). The cell porosity $\varepsilon _f^\ast $ is computed according to the constraint (2.1). Based on $\varepsilon _f^\ast $ and $\varepsilon _{coke}^\ast $, the sub-grid properties are updated, including the permeability ${k^\ast }$, the diffusion coefficient $\hat{D}_{i,eff}^\ast $ and the specific coke surface area $\alpha _v^\ast $;

3. (3) the gas continuity equation is explicitly solved to obtain the new gas density $\rho _f^\ast $;

4. (4) solve the discretized momentum equation (2.6) with the preconditioned biconjugate gradient (PBiCGStab) method to predict the velocity ${\bar{\boldsymbol{U}}^\ast }$ and the mass flux $\textrm{ph}{\textrm{i}^\ast }$. Note that the predicted fields ${\bar{\boldsymbol{U}}^\ast }$ and $\textrm{ph}{\textrm{i}^\ast }$ do not satisfy the mass conservation after this prediction step;

5. (5) the evolution of O₂ mass fraction, (2.9), is solved with the implicit source term scheme to update ${\bar{\omega }_{f,{\textrm{O}_\textrm{2}}}}$. With the new $\bar{\omega }_{f,\textrm{C}{\textrm{O}_\textrm{2}}}^\ast $, the coke reaction rate is recalculated for the following CO₂ equation to obtain a new $\bar{\omega }_{f,\textrm{C}{\textrm{O}_\textrm{2}}}^\ast $. In this step, each discretized equation for the species transfer is solved with the PBiCGStab method;

6. (6) the thermophysical properties are updated based on the new mass fractions of species compositions. The energy equation, (2.15), is implicitly solved with the PBiCGStab method to compute the new temperature ${\bar{T}^\ast }$;

7. (7) update the density $\rho _f^{{\ast}{\ast} }$ and the compressibility ${\beta ^\mathrm{\ast }}$ with the new temperature ${\bar{T}^\ast }$. Then, the pressure equation is implicitly solved with the geometric-algebraic multigrid (GAMG) method to correct the pressure ${\bar{p}^{{\ast}{\ast} }}$ and update the velocity ${\bar{\boldsymbol{U}}^{{\ast}{\ast} }}$ and mass flux $\textrm{ph}{\textrm{i}^{{\ast}{\ast} }}$. After this inner corrector step, the density $\rho _f^{{\ast}{\ast} \ast }$ is also recalculated with the new mass flux $\textrm{ph}{\textrm{i}^{{\ast}{\ast} }}$ through the gas continuity equation, (2.7). This inner corrector step can be looped several times to improve the result accuracy. In this work, 1–2 loops are suggested;

8. (8) check if the PIMPLE algorithm converges at the current time step or reaches the maximum number of outer correctors. The criteria for the time step convergence are defined as absolute tolerances of the energy equation and species equations, all of which should reach 10⁻⁶ in this work. The maximum number of outer correctors is set as 15 with regard to the numerical stability. If the time step convergence or maximum outer corrector number reaches the specified conditions, the time moves onto the next step, otherwise it loops over the entire equation system again.

The adjustable time step size is established by the maximum Courant number, the minimum chemical time scale and the maximum allowed time step of 10⁻⁴ s. Considering the excellent number stability of the PIMPLE algorithm and the fact that the fluid flow in the coked porous medium does not rapidly vary with time, the maximum Courant number can be more significant than 1 and as large as 100–1000 for the present problem. Additionally, the much greater time step over the chemical time step introduces numerical instability, particularly for coke combustion at high temperatures. According to the numerical experiments, the recommended time step is less than 70 times the minimum chemical time step to balance the numerical accuracy, stability and efficiency.

3.3. Validations

The developed micro-continuum models were validated with benchmark problems to demonstrate that the numerical models can accurately simulate complex coke combustion coupled with conjugate heat transfer, heterogeneous reactions and solid morphological evolution. These three principal physicochemical processes were progressively tested using three benchmark cases: thermal flow around heated obstacles with conjugate heat transfer; thermal flow around solid obstacles with heterogeneous reaction; and solid dissolution with structural evolution. The numerical results agree with widely accepted solutions. Notably, the benchmark case of thermal flow around reactive obstacles was compared with the LB models (Xu et al. Reference Xu, Long, Jiang, Zan, Huang, Chen and Shi2018b; Lei et al. Reference Lei, Wang and Luo2021), the present pore-scale micro-continuum model and the commercial COMSOL simulator. The comparison indicated that satisfactory consistency could be achieved for the simple problem without structural evolution and multiscale physics. Refer to Appendix B for more details.

3.4. Model improvements

Coke combustion in an ideal porous medium, where coke deposition occurs on circle-like grains, was revisited to demonstrate the improvements in the developed micro-continuum model compared with previous LB models (Xu et al. Reference Xu, Long, Jiang, Zan, Huang, Chen and Shi2018b; Lei et al. Reference Lei, Wang and Luo2021). The same initial condition, boundary condition and thermophysical properties (constant density) as in the reference (Xu et al. Reference Xu, Long, Jiang, Zan, Huang, Chen and Shi2018b) were employed using the micro-continuum model; however, different numerical results were observed, especially the combustion temperature. Compared with previous studies (Xu et al. Reference Xu, Long, Jiang, Zan, Huang, Chen and Shi2018b; Lei et al. Reference Lei, Wang and Luo2021), the present micro-continuum model, with improvements, can yield reasonable coke combustion phenomena, including two-way coupling between the chemical reaction and mass flows, improved energy conservation during geometrical evolution, and reactive transport in the sub-resolution effective continuum. As a result of the inherited excellent computational efficiency of the micro-continuum framework, the image-based computational domain and the physical duration can be further expanded for more physics, such as reactive fingering. Comprehensive comparisons with the previous model can be found in the supplementary material.

Here, the significance of reactive transfer in unresolved nanopores was discussed to highlight the necessity of image-based simulation of coke combustion through a multiscale porous medium. For reactive flow in the sub-resolution effective continuum, the transport and geometrical properties, particularly the specific reactive surface area, should be carefully considered. The diffuse interface model has been widely applied to numerical simulation of the solid–liquid dissolution problem, which captures the sharp reactive interface between the acid fluid and low-porosity low-permeability low-diffusivity rocks. Considering previous investigations of the coke combustion problem (Xu et al. Reference Xu, Long, Jiang, Zan, Huang, Chen and Shi2018b; Lei et al. Reference Lei, Wang and Luo2021), the heterogeneous coke reaction with O₂ was modelled as the interface condition, which also presumed that the coke surface reaction only occurred at the stationary exterior interface between the coke and fluid cells. However, the significant difference between coke combustion and rock dissolution is that the coke phase is highly porous with numerous nanopores, as shown in figure 1. These nanopores permit O₂ to penetrate primarily through Knudsen diffusion and react on the internal pore surface so that reactive transfer through the sub-resolution pores should be considered. Additionally, the coke pore structure changes significantly during combustion, which may affect the coke combustion characteristics. Thus, it is beneficial to account for the coke pore structure using the random pore model instead of the diffuse interface model for the coke combustion problem. This is commonly accepted in the kinetics modelling of gas–carbonate solid reactions (Bhatia & Perlmutter Reference Bhatia and Perlmutter1980).

By comparison, the coke combustion diffuse surface model and random pore model were simulated to emphasize the difference with/without the sub-resolution reactive surface area and thus the sub-grid chemical reactions. Contours of the porosity, O₂ concentration, reaction heat rate and temperature for these two cases at time instant t = 12.0 s are shown in figure 4. In figure 4(a), the combustion dynamics in the chemical kinetics-influenced regime are shown for the case with the diffuse interface model, where the characteristic chemical time scale was longer than the transport time scale. In this regime, the oxygen penetrated all of the pore space, including the unresolved nanopores inside the coke. The coke combustion reaction spanned the entire computational domain; however, the reaction was strictly restrained to the exterior coke surface regardless of the active internal surface, which had been exposed to the available oxygen. Owing to the unphysical limitation of the diffuse interface model, the sub-grid chemical reaction was neglected, and the reactive surface area was greatly under estimated, which led to repressive reactivity and a low burning temperature. In contrast, the case with the random pore model behaved similar to the typical transport-limited process. Snapshots of the porosity and the combustion heat release rate distribution in figure 4(b) show that a narrow flat combustion front arose and divided the domain into a coke-depleted region and a coke-rich region. The oxygen penetrated the coke-depleted region and then exhausted at the combustion front. Coke was observed to react with O₂ at both the exterior and internal surfaces of the coke, which led to increased coke reactivity and a higher burning temperature.

Figure 4. Coke combustion problem for discussions on model improvements in the reactive surface area model: contour comparisons of the porosity, O₂ concentration, reaction heat rate and temperature around the time instant t = 12.0 s for two cases with diffuse interface model (a) and the random pore model (b), with both simulated at an initial temperature of 673 K and air flux of 2.619 sm³ (m² h)⁻¹. The counter-field contours are coloured with the same contour levels. In the porosity contour plot, the coke region with the intermediate porosity is depicted by the pink colour.

For quantitative comparisons, temporal evolutions of the combustion dynamics were recorded for cases with and without sub-grid chemical reactions, as illustrated in figure 5. From the temporal results, the variations in the coke reaction rate and the combustion temperature were greatly underpredicted at the firing stage without sub-grid chemical reactions (diffuse interface model). Driven by the intensive oxygen diffusion transport owing to the inlet effect, the combustion temperature rapidly increased from the initial temperature of 673 K to approximately 735 K in less than 1 s for the case with sub-grid chemical reactions (random pore model); meanwhile, the peak volume-averaged coke reaction rate exceeded 3.0 kg (m³ s)⁻¹. Accelerated coke combustion with sub-resolution chemical reactions can favour the growth of the combustion front and transform the combustion into the transfer-limited regime. As the combustion front proceeded along the channel, the combustion heat release was gradually balanced with the heat loss from the domain boundaries and, thus, the combustion eventually entered the thermal equilibrium stage. At the thermal equilibrium stage, the coke reaction rates for these two cases were comparable and a burning temperature difference of less than 10 K was maintained. Although the burning temperature difference was not significant, completely distinct combustion dynamics were obtained as a direct consequence of integrating the sub-grid chemical reactions or not. Therefore, the present study emphasized the significance of the sub-grid reactive transfer to simulate non-isothermal reactive flow through a multiscale porous medium. For the coke combustion problem in the present study, the random pore model was employed in the following simulation. Furthermore, combustion regimes from two different reactive surface models were further compared with some experimental measurements in § 4.1 to validate the accuracy of the random pore model.

Figure 5. Coke combustion problem for discussions on model improvements in the reactive surface area model (initial temperature of 673 K and air flux of 2.619 sm³ (m² h)⁻¹): (a) the temporal evolution of volume-averaged coke reaction rate and cumulative conversion; (b) the temporal evolution of the maximum transversely averaged temperature and the transversely averaged O₂ molar concentration at the outlet, for cases with the diffuse interface model and the random pore model.