Modelling of Catalytic Combustion in a Deformable Porous Burner Using a Fluid–Solid Interaction (FSI) Framework

Abstract

The various concepts involved in the mathematical modeling of the fluid–solid interactions (FSIs) of catalytic combustion processes occurring within a porous burner are presented and discussed in this paper. The following aspects of them are addressed: (a) the relevant physical and chemical phenomena appearing at the interface between the gas and the catalytic surface; (b) a comparison of mathematical models; (c) a proposal of a hybrid two/three-field model, (d) an estimation of the interphase transfer coefficients; (e) a discussion of the proper constitutive equations and the closure relations; and (f) a generalization of the Terzaghi concept of stresses. Selected examples of application of the models are then presented and described. Finally, a numerical verification example is presented and discussed to demonstrate the application of the proposed model.

1. Introduction

The science of the catalytic combustion of gases is very complex. It developed as a combination of many practical discoveries, laboratory experiments and as a result of industrial experiments. The complexity of the phenomena grew from some tremendous and intuitive prediction combinations of two components: physical responses and chemical reactions. Even fundamental characteristic parameters such as catalytic activity/inactivity, catalytic internal kinetics, the length, and the temperature of activation are variable and difficult to describe within one unified scheme.

On the other hand, modern industrial systems need compact devices for performing combustion. Therefore, in industrial practice, many different porous burners have been developed and realized. From the literature, it is well know how extremely difficult it is to describe—theoretically and numerically—the reactive flows of gases in porous media, with or without the action of catalysis. There is a lack of advanced models for describing these processes. This lack and these difficulties arise from the fact that single laboratory experiments over a solely catalytic surface cannot be directly repeated in porous media conditions. The knowledge about flows of mixtures in porous media cannot be easily and directly used in the formulation of the fundamental characteristics of combustion in porous media.

The reason for this complexity comes from the fact that the huge internal surface areas of porous media change the flow conditions of gases, both in the domain with or without chemical reactions. In particular, these difficulties come from the fact that the chemical and physical properties of the inner surfaces change during contact with different gases. Therefore, surface phenomena such as slip velocity, surface mobility, thermophoresis, and surface diffusion are different for different gases on the same surface material.

Based on our experience, we recognize that three above-mentioned questions overlap within a one domain of science called fluid–solid interaction (FSI). The three main domains of scientific interest are shown schematically in Figure 1. The practical and scientific knowledge on the indicated domains has been developed very quickly, but separately, over the last two–three decades. Therefore, one can find numerous discrepancies and incoherencies when looking for more adequate mathematical models of catalytic combustion in porous burners. In particular, there is a lack of robust models when the precise mathematical modeling also needs a detailed description of the solid material, such as thermal deformation, sintering, hydration, sorption, etc. Thus, turning our attention to this problem, we supplemented the phrase “porous burner” in the title. Instead of “porous burner”, we wrote “deformable porous burner”.

Hopefully, the common framework for the integration of the three above-mentioned scientific principles into one common mathematical model is given by the concept of fluid–solid interactions (FSIs). FSIs are important in every case when the properties of thin layers of solid–fluid contacts are dominant within the whole phenomena of a model. Continua such as foams only contain a continuum of thin solid–fluid interfaces. Therefore, within the process of modeling continua with a very packed density of the interface surfaces, it is important to average the local surface properties into parameters of a bulk three-dimensional porous medium.

1.1. Catalytic Combustion and Experimental, Theoretical, and Numerical Research

Extensive experimental and theoretical attention has been paid to catalytic combustion over the past three decades. A lot of evidence for the significant potential of heterogeneous processes in reducing the emission pollutants, improving ignition and enhancing the stability of the generated flames has recently been recognized.

Among the many possible examples of experiments there are to take as patterns of modeling these processes, there are two baseline configurations that are often used to experimentally investigate catalytic combustion: stagnation flow fields over a catalytically active foil [1,2,3] and chemical reactors with a catalytically active wire inside them [4,5,6]. In both setups, the temperature of the catalyst is controlled by the resistive heating of the catalytic foil or wire. It is important that the mathematical modeling and simulation of such heterogeneous systems include the coupling applications between the reactive flows and the surface to capture the interactions between them at the interface. Therefore, computational tools for both systems have been recently developed, providing the ability to analyze the elementary chemical and transport processes at the gas–surface boundary and to couple them with a description of the surrounding gas phase. The aim of the undertaken efforts was to achieve a quantitative understanding of catalytic combustion (see Badur et al. [7], Badur and Ochrymiuk [8], Deutschmann [1]).

It should be noted that high embodied energy gases such as hydrogen and high-emmision gases like methane are commonly used in practice. However, recently, carbonless gases have become very important (for instance NH₃). Therefore, from the point of view of the fundamental science, the various interactions of gases with catalysts are currently being studied at the interaction level of the gas phase molecules coming into contact with the elements of the crystal lattice of solid-phase catalysts [9].

Great progress has been achieved in this direction. In particular, phenomena such as bi-stability, surface waves, spatiotemporal chaos, patterns, chemical turbulence, and others have been studied [10].

Catalytic combustion frequently involves some reaction enhancement. An enhancement analysis was developed as an addition to the surface reaction mechanisms in order to find the reaction-controlling steps. According to practical [6,11], catalytic combustors are ceramic honeycomb-shaped devices coated with platinum or palladium. These elements are engineered to maximize the catalytic reaction—a hexagon shape ensures the maximum surface area for the minimum material usage. Strategically, catalytic combustors are designed to burn difficult or incombustible particles that are situated inside the smoke path. Catalytic combustors literally cause “difficult gases” to burn as fuel, creating more heat from the catalytically enhanced reactions. This means that dedicated fuels, such as dangerous gases, can be fully burned, transforming most of their internal energy into useful heat instead of them being released to the environment as pollution.

Yet another enhancement effect can be obtained by using perovskite catalysts with matrix-stabilized combustion in porous ceramic media [12]. If highly porous silicon carbide ceramics are used as the porous media, a catalytically enhanced super-adiabatic combustion of a lean mixture of methane and air can be easily performed. Robayo et al. [12] were also able to provide a direct regulation of the combustion flame and to obtain the effect of “perovskite catalytic enhancement of SiC” due to the use of a specially designed stainless-steel deformable chamber incorporating a quartz window.

However, from the point of view of the environmental impact of combustion, catalysis is used in the most popular catalytic converters—their application in the automotive industry has helped to reduce the emission of pollutants significantly. Special attention needs to be paid to the development of micro-catalytic combustors using high-precision ceramics. It should be noted that a catalytic burner is not a filter. Instead of trapping unburned particles, the combustor deploys chemical catalysis to break them apart. In particular, platinum and palladium atoms loaded into honeycomb-shaped cells trigger their chemical reaction at the contact surface.

A micro-scale catalytic combustor fueled by butane was investigated in the literature [13,14]. In particular, Okamasa et al. [3], developed a high-precision ceramic-tape-casting technology, which was applied in a three-dimensional combustor structure with embedded heat exchange channels. Nano-porous alumina fabricated through the anodic oxidation of aluminum layers was employed as the support for a Pd catalyst. Combustion experiments were carried out in a solder bath to keep the catalyst temperature constant. Complete fuel conversion for an n-butane flow rate of 5.0 sccm was achieved at 390 °C, corresponding to a 100 MW/m³ heat generation. The reaction constants for the catalytic combustion on the Pd/nano-porous alumina were determined with the aid of a 1D plug-flow model. It was also shown in a preliminary experiment in air that the reaction could be self-sustained at 425 °C with an n-butane flow rate of 15 sccm [3].

There have been a number of studies on the catalytic effects and benefits of various materials in terms of enhancing the combustion of hydrocarbon fuels. These studies start from the paper by Karim and Kibrya in 1986 who experimentally investigated the lean blowout limit of the combustion of methane and hydrogen in air [2]. The burner contained a cylindrical combustor that was 150 mm in diameter with a metallic wire mesh acting as a porous media. They coated the wire mesh with eight different materials by using an electroplating method. The materials listed in decreasing order of effectiveness are Pt, Cu, Ag, Brass, Cr, Cd, Ni, and stainless-steel 36, with the Pt coating supporting the lean combustion of methane down to 2.7% by volume, a 0.26 equivalence ratio, and the lean combustion of a 50% methane/hydrogen mixture down to a 0.15 equivalence ratio.

When comparing the obtained results against hydrogen combustion, it was found that hydrogen is more sensitive to the catalytic effects of the materials than methane, particularly at lower temperatures. In order to minimize the thermal stress at high operation temperatures, the catalyst arrangement was redesigned to the reduce temperature gradients. The surface reactions in the micro channels led to stable flames and an extremely high heat generation density of 2–3 ×10⁸ W/m³.

In addition, other authors (Dupont et al. [4]) investigated the stability of methane/air combustion, as well as the emission of pollutants and the radiation efficiency in a honeycomb-shaped porous media containing platinum and palladium catalysts. The honeycomb was made of cordierite and had 400 square cells per square inch. They found that the palladium catalyst supported a lower inlet concentration than the platinum and that it had a proportional effect over all the inlet flow rates tested. The results showed that the minimum stable molar concentration of CH4 was 4.4%, corresponding to 4.4 grams of Pd per piece and a 60 L/min flow rate [4].

1.2. Nickel Catalyst

Palladium is not the only chemical element used as a catalyst. Many companies have based their products on nickel instead (see Figure 2). Nickel-based catalysts exhibit an extremely high catalytic activity in methanol decomposition and in the synthesis of other gasses, and they promote the production of carbon nanostructures (mainly carbon nanotubes) (Badur et al. [7], Jóźwik et al. [5]). One of the most common Ni-based solid-state catalysts is intermetallic-phase Ni₃Al and its alloys, which belong to a family of multifunctional materials, combining the properties of both structural and functional materials. According to [15], intermetallic Ni₃Al thin foils exhibit extremely high catalytic properties in hydrocarbon decomposition reactions.

On the other hand, the relatively high temperature required for maximal hydrocarbon conversion is the main disadvantage of this material. Nevertheless, the high temperatures of the process can be utilized by placing a regenerative heat exchanger downstream of the reactor. Jóźwik et al. [5], proposed a design for a thermo-catalytic reactor with thin strips or foils based on intermetallic-phase Ni₃Al, which appeared as an innovative and extremely promising technology. An example of an alloy foil package based on intermetallic-phase Ni₃Al constructed as a rolled-up sine structure is shown in Figure 2.

1.3. Combustion in Porous Media

It is well known that the conventional burning techniques of fuel–air mixtures have defined flammability limits, beyond which a flame cannot self-propagate due to heat losses. Porous media combustion involves two main scenarios: surface- and matrix-stabilized combustion, which are primarily defined by the flame location in each case. With surface-stabilized combustion, a flame sheet is developed on the surface of a solid porous body by many small individual laminar premixed flames. With matrix-stabilized combustion, the combustion process takes place within the solid porous media. Matrix-stabilized combustion in a porous medium is an advanced technique in which a solid porous matrix within the combustion chamber accumulates heat energy from the hot gaseous products and preheats the incoming reactants. This heat recirculation extends the standard flammability limits and allows the burning of ultra-lean fuel mixtures, conserving energy resources, or the burning of gases with a low calorific value, utilizing otherwise wasted resources. In matrix-stabilized combustion, a low-porosity inlet section is used to transfer heat from the combustion chamber to the reactants and to prevent upstream ignition. The low-porosity inlet section has a pore diameter less than the flame-quenching diameter at the operating conditions, which prevents flashback occurring, which is where the flame speed is higher than the mixture velocity and the flame propagates upstream [16].

The modelling of combustion within porous burners needs to take into account the properties of the material and to modify the momentum balance for the reacting gas mixture, and two energy balances are required for determining the two fields of temperature—the gas mixture and the solid. The optimal permeability for the preheating zone is usually identified following the so-called Carman–Cozney permeability model (Sobieski and Trykosko [17]). The modelling of pore sizes is important if the pores are smaller than the flame-quenching distance in order to prevent the flame from propagating outside of the central section (Krakowska et al. [18]).

Some studies are also available on the addition of a catalyst to a porous material for improving the properties of burners, such as studies on radiant efficiency and the emissions of existing designs. Of these, few have combined these changes with a catalyst material. With regards to materials for porous radiant combustion, both ceramic and metal compositions have been extensively explored, and a review of the most common materials is provided in the paper by Tierney and Harris [19]. For ceramics, the most common materials include cordierite, mullite, alumina, silicon carbide, zirconia or combinations of these. Ceramic materials demonstrate excellent temperature stability, making them attractive for use in porous burners [20].

1.4. Mathematical Modeling of Micro Flows

The main problems that are related to micro flows were described in a paper by Badur et al. [21,22,23]. The description of these phenomena includes the Stefan slip velocity, surface turbulence, surface mobility, surface jumps in temperature, and other parameters. This means that, as has been observed within experiments with a reacting gas and a catalytic surface, it is important to know the detailed relationship between the hydrodynamic and chemical processes on the catalyst surface.

Unfortunately, due to the usually complex composition of catalyst materials, it is difficult to theoretically predict the chemical reaction constants (Arrhenius and so on), and there are no publications concerning numerical research on catalysts based on, for instance, palladium-based coatings. The surface layer depends significantly on the solid temperature, on the strain rate in the near-wall boundary layer, and on the composition and temperature of the reagent flow. It is possible that at relatively low surface temperatures and low hydrogen concentrations in the external flow, heterogeneous reactions can predominate, while at higher temperatures and higher concentrations, homogeneous reactions can predominate. We refer to these phenomena as thermochemical inertia.

One example of measuring thermochemical inertia, which is governed by the time characteristic of the kinetics of surface catalysis, was performed by Kim et al. [9]. This time characteristic of catalytic kinetics on a palladium-based surface is much less than the time characteristics of convective or diffusive processes. Nevertheless, generally, the relationships between micro-flow characteristics, such as the Stefan slip velocity, and micro-catalytic reactions have not yet been determined and need further research [1].

1.5. FSI Framework Description

Regarding the mathematical modeling of catalytic combustion, the equations must take several phenomena into account in order to rigorously simulate a single porous channel. The key points to consider include:

(1)

Heterogeneous reactions at the catalyst surface layer and homogeneous reactions in the gas phase;

(2)

Mass, momentum, and thermal energy transfer by convection and diffusion in the gas phase and at the gas–solid catalytic layer;

(3)

Axial heat transfer in the solid phase through the active catalyst and substrate by conduction and radiation. This assumption requires the modeling of three fields of temperature, namely the gas, solid, and catalytic layer;

(4)

In addition, these phenomena are interact strongly because of the intense thermal effects associated with the heat of combustion released within the catalytic layer. These require a simultaneous simulation of porous solid thermal deformation that changes the distribution of the porosity.

The available literature has approached the phenomena of catalytic combustion from many angles (see a review by Deutschmann [1]). In the authors’ group, the so-called FSI framework has been developed for many years (Badur et al. [21,24]; Ziółkowski and Badur, [25]; Ziółkowski et al. [26]), especially for porous materials (Badur et al. [22], Ochrymiuk [27]) and thermal–FSI coupling (Froissart et al. [28,29]; Karcz and Badur [30]; Kraszewski [31], Ochrymiuk [32], Ziółkowski and Badur [33]). The great advantage of this approach is the simultaneous description of both of the unknown fields, which reduces the necessity of iteratively coupling the physical quantities between the fluid and solid.

1.6. Main Goals of the Article

The starting point of the presented paper is the model proposed by the authors’ group in previous articles (Badur et al. [7], Badur and Ochrymiuk [8]), which were focused on the heat transfer within a micro catalytic combustion process. This is based on a detailed chemical reaction mechanism, including a gas-phase reaction and a surface catalytic reaction. As an extension to this, a new multi-field model is proposed for adopting FSI software to simulate a real catalytic combustion process of a premixed hydrogen–oxygen flow in a sub-millimeter ceramic/alloy porous burner.

The second section describes the subject of continuum mixture formation, both in the fluid and the solid. The third section is a crucial one, because it deals with the description of the behavior of the catalytic layer in the contact zone between the fluid and solid. The equations of mass, momentum, and energy conservation within the boundary of this catalytic layer are presented. The novelty of this model is driven by the detailed modeling of the surface geometry and surface kinematics.

In Section 4, an averaged two/three-field equation for a porous catalytic combustor is presented. The concept is based on the introduction of the volumetric (not surface, so far) coefficients of mass, momentum, and energy exchange.

Section 5 contains the steps taken in developing the kinematic relationships involvedp. Section 6 describes the volumetric transport closure. Section 7 presents various new constitutive relations. Section 8, Section 9 and Section 10 are devoted to the presentation of the advanced concept of Terzaghi stresses and the concept of the material effort involved in ceramic composites.

In Section 11, the concept of taking numerical readings of disaggregated energy production is shortly presented—it has some similarity with concepts of Bejan [34], and Feidt [35,36]. Finally, the last section contains an example of an analysis based on the concept of thermal FSI.

2. Governing Equations

2.1. Evolution of Species within the Gaseous Mixture

The set of governing equations involved in this process contains the evolution of $N_s-1$ species, where the production rate of a species increases or decreases subject to the chemical reaction below (Badur and Ochrymiuk [8], Truesdell [37]):

$${\partial }_t\left(\epsilon \rho Y_i\right)+\mathrm{div}\left(\epsilon \rho Y_i{v}_g\right)=\mathrm{div}\left(J_i\right)+\epsilon {\dot{\omega }}_i{W}_i, i=O_2,N_2,C{H}_4, H_2,\dots ,N_S-1$$

(1)

where the diffusive flux of the i-th species is denoted by $J_i$, which is defined by the Dixon–Levis law:

$$J_i=\epsilon \rho D_{mi} \mathrm{grad} \left(Y_i\right)$$

(2)

where $W_i$ is the molar density of each species. The rate of production of every species, ${\dot{\omega }}_i$, which is related to a single mole of the mixture, depends on the intensity of the $N_R$ chemical reactions, the stoichiometric coefficients $v_{ki}^{\prime }$, $v_{ki}^{″}$, and the chemical progress coefficients $k_k^f, k_k^r$, which depend on the Arrhenius expression (Badur and Ochrymiuk [8]):

$$k_k^t=A_k{T}_g{}^{{\beta }_k}\exp \left(\frac{-E_k}{R T_g}\right), t=f,r$$

(3)

The activation energy, $E_k$, coefficients $A_k, {\beta }_k$, and the stoichiometric coefficients were obtained from the chemical tables in GRI-71. The gas constant of the mixture is the sum of the species: $R={{\displaystyle \sum }}_{i=1}^{N_S}{Y}_i{R}_i$.

2.2. Evolution of Species within the Composite Porous Material

Ceramic matrix composites (CMCs) are a subgroup of composite materials and a subgroup of ceramics. They consist of ceramic fibers embedded in a ceramic matrix. Both the fibers and the matrix can be manufactured from any ceramic material, including carbon and carbon fibers, which also can be regarded as a ceramic material. Recent efforts to build a better mathematical composite model have been supported by some similarities with models of mixtures of reacting gases. The first element is familiar only in terms of the physical properties of the constituents (matrix and fibers), and the sintering process is analogous to the reaction processes in gases. It is evident that a simple procedure of adding together the properties of the mixture is wrong from the very beginning, thus this analogy improves the predictions diametrically. Therefore, in the modern literature, researchers have focused on obtaining proper geometrical data before modeling the sintering process. These data can be obtained from image analysis processes obtained via scanning electron microscope (SEM) scans of the geometrical and compositional properties of the tested materials. The volume fraction of the fibers is chosen during the preparation of the material and usually makes up 1–5% of the whole volume. Regardless of the chosen experimental method, the input data for meso-mechanical analysis are not always straightforward to measure. Multi-scale strategies are helpful in generating elegant homogenized mesoscopic properties, which can be numerically identified from simulations performed at the micro-scale (Ochrymiuk [32], Orłowski et al. [38,39]).

In addition, some homogenization procedures are needed to take into account the sintering of two ceramic components. This step is much more difficult than the procedure of homogenizing heterogeneous reactions (see Equation (1)) of gases. Usually, the Mori–Tanaka model of sintering and double homogenization models for representative volumes are utilized within professional software in order to account the different physical characteristics of the composite material. Unfortunately, another factor needs to be taken into account, as compound ceramics have a distinct weak interface phase or a layered phase that is distributed throughout the bulk ceramic matrix. By varying the type and distribution of the weak interface phase or the layered phase in the composite it is possible to obtain a wide range of mechanical properties [40]. Due to limited space, this problem is not discussed in this paper; however, it must be solved within the same FSI frame.

3. FSI—The Governing Equations of Fluid–Solid Contact

3.1. The Concept of a Catalytic Surface Layer

Catalytic combustion needs the presence of a fuel–oxygen mixture contacting a the solid surface that is covered, for instance, with palladium. As experiments have shown, the combustion of hydrogen/oxygen mixtures on palladium exhibits bi-stability for lean mixtures. This happens due to the presence of a special thin gas layer that holds special physical–chemical properties, such as Stefan slip velocity [1,11,41]. Therefore, the concept of a thin catalytic interface layer formed with a special gas is correct, and so this can be used to apply the proper boundary conditions of a fluid–solid contact surface [42]. This is especially important when dealing with porous solid composite materials, where the area of the contact surface is extremely high. This paper intended to expand the range of mathematical modeling methods based on the framework of the FSI (fluid–solid interaction) approach, which describes complex continua. The basis for this further development was papers published within the Energy Conversion Department of the Institute of Fluid-Flow Machinery at the Polish Academy of Sciences (Badur et al. [23], Ziółkowski et al. [43], Ziółkowski and Badur [25], Ziółkowski et al. [26]).

The list of the main parameters characterizing catalytic surface layers includes: temperature $T_{cat}$, density ${\rho }_{cat}$, special surface velocity (called the Stefan slip velocity, $s$), and the coverage by the absorbed species ${\phi }_{\left(b\right)}$. Reactive flow with catalytic properties needs to be coupled with the heterogeneous chemical reactions and the transport equations at the gas–surface interface. Therefore, masses and energy conservation equations for the species need to be established at the interface, considering a small control volume with a finite thickness of the catalytic layer.

The geometry of a catalytic layer is shown in Figure 3. It can be characterized by the gas layer density ${\rho }_{cat}$ (kg m⁻²), the particle layer velocity $s$ (m s⁻¹), and the surface momentum ${\rho }_{cat}s$ (kg m⁻¹ s⁻¹). In addition, the surface excess of the flux of momentum $t_{cat}$ can be incorporated [21].

The transport of momentum within the surface layer can be slightly different from the bulk of a typical gas mixture. For instance, it can contain recoverable (elastic) transverse components. Then, if $a_{\alpha }, n$ ($\alpha =1,2$) are the base vectors on the middle surface of the layer $\mathsf{\Sigma }$, where $n$ is the unit normally oriented to the surface $\mathsf{\Sigma }$, then the surface momentum flux has complicated components:

$$t_{cat}=t^{\alpha \beta }{a}_{\alpha }⨂a_{\beta }+t^{\alpha n}{a}_{\alpha }⨂n+t^{n\beta }n⨂a_{\beta }+t^{nn}n⨂n$$

(4)

This is also called the surface Cauchy stress tensor. The physical properties of the layer are unknown a priori since they depend on the resulting apparent properties of both of the contacting continua. For instance, most simple situation is the case of water contacting with air. In this case, $t_{cat}$ reduces to the surface tension, which is a two-dimensional spherical tensor tangentially oriented to the surface: $t_{cat}=\gamma I_2$, where $\gamma$ is the water–air surface tension [kg/msm] and $I_2=I-n⨂n=a^{\alpha \beta }{a}_{\alpha }⨂a_{\beta }$ is the surface metric tensor [21,44,45]. A second metric surface tensor is usually called the curvature tensor and can be defined as the surface gradient of a normally oriented vector:

$$I{I}_2=-{\mathrm{grad}}_2\left(n\right)=-n⨂{\nabla }_2=b^{\alpha \beta }{a}_{\alpha }⨂a_{\beta }$$

(5)

This has the following well known variables: $I_b=\mathrm{tr}\left(I{I}_2\right)=b_1^1+b_2^2=r_1^{-1}+r_2^{-1}$, which is the mean curvature, and $I{I}_b=\det \left(I{I}_2\right)=\det \left(b_{\alpha \beta }\right)=r_1^{-1}{r}_2^{-1}$, which is the gaussian curvature. For example, $I_b$ appears in the seminous Laplace formulae.

The gradient of the surface vector, for instance, the gradient of the Stefan slip velocity $s=s_{\alpha }{a}^{\alpha }+s_{n}n$, can be calculated as follows:

$$\begin{array}{c}{\mathrm{grad}}_2\left(s\right)=\left(s_{\alpha }{a}^{\alpha }+s_{n}n\right)⨂{\nabla }_{\beta }{a}^{\beta }\\ =\left(s_{\alpha ;\beta }-s_n{b}_{\alpha \beta }\right)a^{\alpha }⨂a^{\beta }+\left(s^{\alpha }{b}_{\alpha \beta }+s_{n,\beta }\right)n⨂a^{\beta }\end{array}$$

(6)

Having obtained the surface gradient, the surface divergence can be calculated as:

$${\mathrm{div}}_2\left(s\right)\equiv C_{1,2}\left({\mathrm{grad}}_{2}s\right)=\left(s_{\alpha ;\beta }-s_n{b}_{\alpha \beta }\right)a^{\alpha \beta }={\mathrm{div}}_2\left(I_{2}s\right)-I_b{s}_n$$

(7)

Analogically, the surface divergence [${\mathrm{div}}_2 t_{cat}]$ can be calculated in the surface momentum balance. Considering the viscous flow within the surface layer, the proper rate of the surface deformation tensor must be defined as a symmetrical part of the surface gradient of the velocity, i.e., $d_2=0.5\left({\mathrm{grad}}_2\left(s\right)+{ \mathrm{grad}}_2\left(s\right)\right)$. Then, the surface viscous stresses, analogously to 3D coordinates, $\tau =2\mu d$, can be written as ${\tau }_2=2{\mu }_2{d}_s$, where ${\mu }_2$ is the Navier surface viscosity coefficient [46,47]. As indicated in the literature, the Navier viscosity is a physically complex concept, since the anisotropy in the surface viscous stresses is usually experimentally observed [48,49]. This leads to the conversion of the single isotropic Navier coefficient, ${\mu }_2$, into a full constitutive expression for anisotropy. For instance, ${\tau }_{2\alpha \beta }={\mu }_{2\alpha \beta }^{ \gamma \delta }{d}_{\gamma \delta }$ ($\alpha ,\beta ,\gamma ,\delta =1,2$), where ${\mu }_{2\alpha \beta \gamma \delta }$ is an anisotropic surface viscosity coefficient tensor [50,51,52].

3.2. Surface Balance of Mass

Generally, the catalytic surface density is a function of the concentration of the catalytic surface species, $Y_{\left(b\right)}$, which is governed by the developing of the species equation, i.e., Equation (18). Nevertheless, similar to the bulk stream, the total mass balance of the catalytic layer can be formulated as:

$${\partial }_t\left({\epsilon }_s{\rho }_{cat}\right)+{\mathrm{div}}_2\left({\epsilon }_s{\rho }_{cat}s\right)={\displaystyle \sum }_{\left(b\right)}^m{\dot{m}}_{gs\left(b\right)}\left(Y_{g\left(b\right)}-Y_{s\left(b\right)}\right)$$

(8)

The sources of mass in Equation (8) cannot be taken arbitrarily; they must fulfil the mass flux fouling equation on the catalytic surface:

$${\displaystyle \sum }_{\left(b\right)}^m{\dot{m}}_{gs\left(b\right)}\left(Y_{g\left(b\right)}-Y_{s\left(b\right)}\right)={\displaystyle \sum }_{i=0}^{N_g}{\dot{\omega }}_i{Y}_i$$

(9)

3.3. Surface Momentum Balance

The starting point for the momentum balance within the catalytic layer is quite a fundamental contact boundary condition for fluid and solid stresses:

$$t_g{n}_g+{\sigma }_s{n}_s=0$$

(10)

where the Cauchy stress in the fluid is $t_g=-p_{g}I+{\tau }_g+{\sigma }_g$ and, similarly, the Cauchy stress in the solid is ${\sigma }_s={\sigma }^{\prime }+\alpha p_{g}I$ (see Equations (48) and (50)). Note that $n=-n_s=n_g$.

This boundary equation can be extended by adding the surface resistivity force (Coulomb, (1801) [53]), $f_r=f_r^{Du}+f_r^{Na}+f_r^{Bu}$, which consists of three parts: the Duhem force (1898) [54], the Navier force (1822) [46], and the du Buat force (1796) [33]. These are determined via the surface viscosity coefficients. The dimensionless forms of these coefficients are called the Duhem number, the Navier number, and the du Buat number (Ziółkowski and Badur [33]). In the case of porous media, where the contact surface area is enormous, the Navier number is more important than the Reynolds number (dimensionless bulk viscosity) (Reynolds, 1901 [55]).

Furthermore, Equation (6) can be extended by incorporating the so-called surface mobility forces, $f_m=f_m^{Ga}+f_m^{Re}+f_m^{Ma}+\dots$, that consist of several parts: the Graham surface separation force (1848) [33], the Reynolds thermal transpiration force (Reynolds, 1879 [56]), the Maxwell transpiration force (Maxwell, 1879 [57]), etc., (Badur et al. [22], Ziółkowski and Badur [25]).

Then, adding the forces within the surface control domain together, $dV=l\times dA$ (Figure 3), the following equation can be obtained:

$$\left(t_g-{\sigma }_s\right)n+f_r+f_m=0$$

(11)

The static balance of the surface forces can be complemented with the surface divergence of the so-called Laplace tensor of momentum, $t_{cat}$ (see Equation (4)), so the static balance of the catalytic layer is governed by:

$${\mathrm{div}}_2\left(t_{cat}\right)+\left(t_g-{\sigma }_s\right)n+f_r+f_m=0$$

(12)

Finally, Badur et al. [21,24], proposed a description of the whole dynamics of the surface layer with the following balance:

$$\begin{gathered} {\partial }_t\left({\rho }_{cat}s\right)+{\mathrm{div}}_2\left({\rho }_{cat}s⨂s\right)+{\mathrm{div}}_2\left(t_{cat}\right)+ \\ +\left(t_g-{\sigma }_s\right)n+f_r+f_m+{\dot{m}}_g\left(v_g-s\right)=0 \end{gathered}$$

(13)

Using the reasoning of d’Alembert and Euler, the definition of a surface acceleration vector is as follows (D’Alembert [58]):

$$a_2=\frac{d_2}{dt}s={\partial }_{t}s+\left({\mathrm{grad}}_{2}s\right)(I_{2}s)$$

(14)

The momentum balance of the layer (Equation (10)) becomes a nonlinear differential equation for two additional unknown fields: the surface mass density ${\rho }_{cat}$ and the layer slip velocity. If $t_{cat}=\gamma I_2$, $t_g=p_{g}I$, ${\sigma }_s=p_{s}I$, then the layer momentum (Equation (9)) leads to the generalized Young–Laplace equation [47,59,60]:

$${\mathrm{div}}_2\left(\gamma I_2\right)+\left(p_g-p_s\right)n=0$$

(15)

In general, the flux of the layer momentum, $t_{cat}$ (see Equation (13)), is responsible for the recoverable (elastic) and viscous transport of the surface momentum, i.e., $t_{cat}=t_{cat}^{\left(e\right)}+t_{cat}^{\left(v\right)}$. The first and most important part of elastic flux is known as the capillarity diade, $\gamma I_2$. However, Stokes, in 1845, introduced an additional “normal” part, $\varpi$, and in 1876, Gibbs added the curvature tension part, $C$, such that [23,44]:

$$t_{cat}^{\left(e\right)}=\gamma I_2+CI{I}_2+\varpi n⨂n$$

(16)

The viscous properties of the catalytic layer depend on the so-called “apparent viscosity”, which, in general, possesses a transversal anisotropy. Using the diade of the surface deformation rate, $d_2=\frac{1}{2}\left({\mathrm{grad}}_2\left(s\right)+{ \mathrm{grad}}_2\left(s\right)\right)$, the following viscous stress tensor can be obtained [25]:

$$t_{cat}^{\left(v\right)}={\lambda }^{\prime }\left(tr d_2\right)I_2+{\lambda }^{″}\left(s_{n,n}\right)n⨂n+2{\mu }^{\prime }{I}_2{d}_2{I}_2+2{\mu }^{″}\left(d_2-I_2{d}_2{I}_2\right)$$

(17)

where the formulas for the four apparent surface viscosities, ${\lambda }^{\prime }, {\lambda }^{″}, {\mu }^{\prime }, and \mu ″$, need further investigation. Nevertheless, these coefficients should be distinguished from the surface viscosity such as the Navier surface viscosity. See the papers by Arkilic et al. [42], Kowalewski et al. [61], Lockerby et al. [52] and Morini et al. [59], where fundamental experiments were developed.

3.4. Catalytic Evolution of Species

The conditions of the catalytic surface are described by the temperature and the coverage with the adsorbed species. The reactive flow has to be coupled with the heterogeneous chemical reactions and the transport at the gas–surface interface. Therefore, conservation equations for species masses and energy can be established at the interface, considering a small control volume of the catalytic layer, $dV=l\times dA$ (Figure 3, adjacent to the surface. Then, the mass fraction of a gas-phase species $\left(b\right)$ at the surface is determined by the diffusive and convective processes and the production or depletion rate $\left(b\right)$ of a given species by surface reactions:

$$\begin{gathered} {\partial }_t\left({\epsilon }_s{\rho }_{cat}{Y}_{\left(b\right)}\right)+{\mathrm{div}}_2\left({\epsilon }_s{\rho }_{cat}{Y}_{\left(b\right)}s\right)= \\ ={\mathrm{div}}_2\left(j_{\left(b\right)}\right)+\left[J_{\left(b\right)}+{\rho }_{cat}s{Y}_{\left(b\right)}\right]\bullet n+{\dot{m}}_{gs\left(b\right)}\left(Y_{g\left(b\right)}-Y_{s\left(b\right)}\right)+{\epsilon }_s{\dot{\omega }}_{\left(b\right)}{W}_{\left(b\right)} \\ \mathrm{for} \left(b\right)=1,2,3,\dots \end{gathered}$$

(18)

where ${\rho }_{cat}$ is the surface density [kg m^-2] of the catalytic layer, $Y_{\left(b\right)}$ is the mass fraction of a species in the control layer, $j_{\left(b\right)}$ is the surface diffusive flux (including thermal diffusion), $J_{\left(b\right)}=J_i$ for $\left(b\right)\equiv i$ species is the mass flux diffusion of species $i$, as discussed in Equation (1). The surface porosity ${\epsilon }_s$ is defined as the ratio of the heterogeneous surface reaction area $A_{cat}$ (desorption minus adsorption) to the unit surface area $A_{geo}$:

$${\epsilon }_s=A_{cat}/A_{geo}$$

(19)

The area, $A_{cat}$, refers to the actual catalytically active surface area and can be determined experimentally, e.g., by chemisorption measurements [14]. The ratio ${\epsilon }_s$ is also used to describe the dependence of the overall reaction rate on catalyst loading and the effects of hydro-thermal aging for structure-insensitive catalysts [62]. It was recently applied to model the performance of on-road-aged three-way catalysts [63].

The surface porosity is very similar to the so-called surface coverage factor ${\phi }_{\left(b\right)}$ i.e., the fraction of the surface sites covered by species $Y_{\left(b\right)}$. The rate of change ${\phi }_{\left(b\right)}$ depends on ${\dot{\omega }}_{\left(b\right)}$ in a simple manner, as described in the paper by Deutschmann [1]:

$${\partial }_t{\phi }_{\left(b\right)}=\beta {\dot{\omega }}_{\left(b\right)} \mathrm{for} \left(b\right)=1, 2, 3,\dots m$$

(20)

where $\beta$ is the surface site density of the catalyst. The numerical solution of the catalytic evolution of a species (Section 3.4) must be performed. As the surface momentum balance (Section 3.3) and the catalytic evolution of species (Section 3.4) are fully time-dependent, transient phenomena such as ignition, oscillation, extinction, or catalyst poisoning can be described in detail. The surface coverage also fulfills the site fraction conservation law, ${\displaystyle \sum }_{\left(b\right)=1}^m{\phi }_{\left(b\right)}=1$.

Returning to Equation (18), the term ${\dot{\omega }}_{\left(b\right)}{W}_{\left(b\right)}$ indicates the source of the heterogeneous surface reaction rate (desorption minus adsorption), which is given in kg m⁻² s⁻¹. The molar net production rate of gas phase species $b$ given as ${\dot{\omega }}_{\left(b\right)}$ (mol m⁻² s⁻¹) now refers to $W_{\left(b\right)}$, which is the actual molecular mass of species $b$; $m$ is the number of gas-phase species on the surface, and $dA$ is the surface layer area. Additionally, in Section 3.4, the outward-pointing unit vector $n$ is normal to the surface, so if chemical surface reactions occur, the slip velocity vector $s$ can be non-zero at the catalytic surface.

This so-called Stefan slip velocity, $s$, which is given by the temperature of the catalyst, is derived from various contributors to the energy balance at the catalyst layer. A non-zero Stefan velocity, $s$, occurs for the net mass flux between the surface and the gas phase [1,64]:

$$s·n={{\displaystyle \sum }}_{\left(b\right)=1}^m{\epsilon }_s{\dot{\omega }}_{\left(b\right)}{W}_{\left(b\right)}$$

(21)

At steady-state conditions, this mass slip disappears unless mass is either deposited on the surface (e.g., chemical vapor deposition) or ablated (e.g., material etching). Equation (6) basically means that for $s=0$, the amount of gas-phase molecules of species $b$, which are consumed/produced at the catalyst by adsorption/desorption, have to diffuse to/from the catalytic wall (Equation (7)). Only for fast transient adsorption/desorption processes, e.g., during the ignition of catalytic oxidation, does the steady-state Equation (6) break down and special treatment of the coupling is needed. Furthermore, these fast transient processes may lead to heat accumulation terms as well as to additional convective transport and the associated pressure gradients in the fluid phase above the catalyst [62].

Concerning ${\dot{\omega }}_{\left(b\right)}$, it is known that the rate of a catalytic reaction is very specific to the catalyst formulation; therefore, only global rate expressions have been used for many years [1,3]. These reaction rates are based on the catalyst mass, catalyst volume, reactor volume, and the catalyst external surface area. The implementation of this surface-kinetics-based approach into mathematical simulations is straightforward; the reaction rate can, in general, be expressed by any arbitrary function of gas-phase concentrations, and the temperature at the catalyst surface, $T_{cat}$, is calculated at every computational cell containing either catalytically active particles or walls.

It is evident that this catalytic-layer-like approach cannot account for the complex variety of phenomena of catalysis and that the rate parameters must be evaluated. The surface chemistry is modeled by elementary reactions on a molecular level in the gas phase and in the solid molecules at the surface. As usual, the temperature dependence of the reaction rate can be described by a modified Arrhenius equation. Special care must be taken for the additional coverage dependence of the rate of some surface reactions (Badur et al. [7]). For reversible reactions, the rate coefficients are related to the forward rate coefficients through the equilibrium constant. For instance, the reactions of hydrogen, oxygen, and methane on polycrystalline platinum along with their rate expressions are given explicitly in the paper by Kang et al. [63]. All the pre-exponential factors, $A$, were chosen to be independent of the temperature. Details on the reaction steps and the rate data are discussed elsewhere (Badur and Ochrymiuk [8], Weisz [14]). The thermochemical data needed to calculate the equilibrium constants for the reversible reactions were taken from [65].

It should be highlighted that the correct and direct mathematical modeling of surface catalytic reactions has up to now been too difficult to be conducted. It is obvious that the direct computation of the chemistry of surface reaction rates at the molecular level leads to a closer comprehensive description, at least for idealized systems. However, powerful methods of mathematical modeling such as direct numerical simulations (DNS), large eddy simulations (LES), lattice Boltzmann models (LBM), density functional theory (DFT), molecular dynamics (MD), and the Monte Carlo (MC) method are currently only in the initial and development stages. Therefore, these approaches cannot be implemented in a complex simulation of flow/deformation fields undergoing a catalytic reaction of technically relevant systems. This is due to a lack of efficient algorithms that are able to reduce the immense amount of computational time currently needed (Deutschmann [1]).

In the presented FSI framework approach, the catalytic system is treated as a black box. There is no alternative; the knowledge gained from experimental and theoretical surface science studies can be implemented in the chemical models used in reaction engineering simulations. A tractable approach involves the treatment of surface chemistry by rate equations that are strongly governed by the surface catalytic temperature $T_{cat}$ and a set of surface coverage factors, ${\phi }_{\left(b\right)}$ (see Equation (6)). They depend on the time and the meso/macroscopic reaction position, but they are averaged over microscopic local fluctuations. The surface structure of the catalyst is associated with a surface site density, $\beta$, that describes the maximum number of species adsorbing on a unit surface area, given in mol m^-2. Each surface species, $Y_{\left(b\right)}$, is associated with a coverage, ${\phi }_{\left(b\right)}$. Under the assumptions made, a multi-step (quasi-elementary) reaction mechanism can be set up. The local chemical molar source term is then defined by the following equation (Badur and Ochrymiuk [8]):

$${\dot{\omega }}_{\left(b\right)}={{\displaystyle \sum }}_{p=1}^{N_s}{\nu }_{\left(b\right)p}{k}_p{\displaystyle \prod }_{\left(j\right)=1}^{N_s+N_r}{\left(X_{\left(j\right)}\right)}^{{\nu }_{\left(j\right)p}^{\prime }}$$

(22)

where $N_s$ is the number of surface reactions, $\left(X_{\left(j\right)}\right)$ is the molar surface concentration (mol m⁻²), and ${\nu }_{\left(b\right)p}$ and ${\nu }_{\left(j\right)p}^{\prime }$ are the stoichiometric coefficients. The expression for the rate coefficient, $k_p$, is still in the research stage, although in the literature there are numerous chemical hypotheses concerning it. The most popular is given by Equation (13) (Badur and Ochrymiuk [8]):

$$k_p=A_p{\left(T_{cat}\right)}^{{\beta }_p} \exp \left[\frac{E_p}{R{T}_{cat}}\right] {\displaystyle \prod }_{\left(b\right)=1}^m{\left({\phi }_{\left(b\right)}\right)}^{{\xi }_{\left(b\right)}}\exp \left[\frac{e_{\left(b\right)}{\phi }_{\left(b\right)}}{R{T}_{cat}}\right]$$

(23)

where two additional coverage dependence factors, ${\xi }_{\left(b\right)}$ and $e_{\left(b\right)}$, appear (Prasad et al. [13]). Here, $A_p$ is the pre-exponential factor, ${\beta }_p$ is the temperature exponent, and $E_p$ is the activation energy. It should be remembered that thermodynamic data for adsorbed species are difficult to measure, so the thermodynamic consistency of surface kinetics for reaction networks is a crucial issue. Two alternative methods have been proposed to enforce thermodynamic consistency without the necessity of explicitly knowing all the thermodynamic adsorbents [13,41].

3.5. Catalyst Temperature—Surface Energy Balance

The temperature of a catalytic reaction, $T_{cat}$, can be determined from a surface energy balance. The key assumption taken forward is that the catalytic layer is able to perform conductive, convective, and diffusive energy transport from the gas phase adjacent to the surface. This includes the chemical heat released at the surface, the thermal radiation, and the resistive heating at the catalyst surface. This assumptions results in the following surface energy equation [33]:

$$\begin{gathered} {\partial }_t\left({\rho }_{cat}{c}_{cat}{T}_{cat}\right)+{\mathrm{div}}_2\left({\rho }_{cat}{c}_{cat}s\right)-T_{cat}{I}_{2}sn+{\mathrm{div}}_2\left(j_{cat}\right)+h\left(T_g-T_s\right)+ \\ +\left(q_g-q_s\right)·n+{{\displaystyle \sum }}_{i=1}^{N_s}{h}_i\left(J_i+{\rho }_g{Y}_{i}s\right)·n+{\sigma }_{em}{\epsilon }_{em}\left(T_{cat}^4-T_{ref,c}^4\right)+ \\ +{\displaystyle \sum }_{\left(b\right)=1}^m{\dot{\omega }}_{\left(b\right)}{W}_{\left(b\right)}{h}_{\left(b\right)}+f·s+I^{2}R=0 \end{gathered}$$

(24)

where $s$ is the Stefan catalytic slip velocity, $f=f_r+f_m c_{cat}$ denotes the specific heat capacity of the catalyst, and ${\rho }_{cat}$ is the density of the catalyst layer. The surface heat flux, $j_{cat}$, is given by an analogy to the Fourier conductivity law, $j_{cat}={\lambda }_{cat}{\mathrm{grad}}_2(T_{cat})$, and $h_{\left(b\right)}$ is the specific enthalpy of the catalytic species $b$, either in the gas phase or at the surface. In the radiation term, ${\sigma }_{em}$ is the Stefan–Boltzmann constant, ${\epsilon }_{em}$ the temperature-dependent surface emissivity, and $T_{ref,c}$ is the reference temperature at which the surface radiates. The term $f·s$ represents the contribution of specific surface work [66].

Additionally, the term $I^{2}R$ represents an energy source corresponding to the resistive heating of the catalyst, where $I$ is the current and $R$ is the electrical resistance, depending on the temperature. In summary, $m$ is the number of surface species, $c_p$ is the specific heat capacity of the gas at the wall, and $h$ denotes a heat transfer coefficient that includes the catalyst volume and a small control volume in the gas phase adjacent to the surface.

4. A Set of the Two-Component Continuum Governing Equations

Once the independent balance equations for gases and solids within the bulk and on the moving boundary (so-called FSI boundary conditions) have been obtained, an efficient model of a porous, two-component continuum where the gas and solid surface take an equivalent role can be constructed. It is obvious that the surface boundary conditions, where divergence [${\mathrm{div}}_2\left(·\right)]$ and surface gradients [${\mathrm{grad}}_2\left(·\right)$] (see Equations (5) and (7)) appear, cannot be further useable in a three-dimensional “equivalent” continuum. Therefore, the dense contributions of those ${\mathrm{div}}_2\left(·\right)$ and ${\mathrm{grad}}_2\left(·\right)$ terms multiplied by the surface density $a$ should be replaced by some “internal closures” that can be interpreted as some internal sources that come from the internal substructure of the continuum mixture.

In order to consider the catalytic combustion process within a porous burner made of modern composites (assuming that the catalytic effect of the porous material at high temperatures is not negligible), any simplification within the porosity effect is a rather wrong assumption, since the “heart” of porous burners are catalytic reactions. The geometry of porous composites can be determined by the following parameters [18]:

-

$\epsilon$—porosity (m³/m³), typical values between 0.7 and 0.9;

-

$a-$ pore surface density (m²/m³), typical values between 1 and 200;

-

$\tau$—tortuosity (m/m), typical values between 1.1 and 2.0.

Please note that in the literature on porous materials [67,68], two-field models dominate.

In the calculation domain, two coordinates dominate, namely $x and r$, or axial and radial, as the circumferential direction in the first approach can be omitted. The porosity of a burner, $\epsilon$, typically has values of $\epsilon =0.7-0.9$. $T_g$ denotes the temperature field of the gases and $T_s$ denotes the temperature of the solid material. In this paper, the well-established European denotations “grad” and “div” are used instead of “nabla”.

The baseline assumption is that the porous mixture possesses a two-component (see Figure 4) type of modelling, namely gas $\left(g\right)$ and solid $\left(s\right)$. As a result, a set of mass, momentum, and energy balances is dedicated to each component.

4.1. Mass Balance

Within the framework of two-field density, two-field momentum, and three-field temperature, there is no more room in the catalytic layer Equation (8). Therefore, in this approach, the transport of mass is governed by a volumetric (bulk) coefficient, ${\dot{m}}_{gs\left(b\right)}$.

4.1.1. Total Mass Balance

The mass balance of the gaseous mixture takes a similar form to the continuity equation of the reacting mixture:

$${\partial }_t\left(\epsilon {\varrho }_g\right)+\mathrm{div}\left(\epsilon {\rho }_g{v}_g\right)={{\displaystyle \sum }}_{\left(b\right)=1}^M{\dot{m}}_{gs\left(b\right)}\left(Y_{g\left(b\right)}-Y_{s\left(b\right)}\right)$$

(25)

the porous material under deformation is given as:

$${\partial }_t\left(\left(1-\epsilon \right){\varrho }_s\right)+\mathrm{div}\left(\left(1-\epsilon \right){\rho }_s{v}_s\right)=-{{\displaystyle \sum }}_{\left(b\right)=1}^M{\dot{m}}_{gs\left(b\right)}\left(Y_{g\left(b\right)}-Y_{s\left(b\right)}\right)$$

(26)

where the velocity of the gas mixture is $v_g$, and the velocity of the solid is $v_g={\partial }_{t}u=\dot{u}$, where $u$ denotes the displacement field of the solid composite. The explicit formula for the density of the gas mixture, ${\rho }_g$, uses the mass, $Y_i$, or volume fraction of species $X_i$ for expressing ${\rho }_g={{\displaystyle \sum }}_{i=1}^{N_S}{Y}_i{\rho }_i$. The index $i=O_2,N_2,C{H}_4, H_2,\dots ,N_S$ deals with the mixture of the species, and $N_S$ indicates how many species there are.

The total transport of mass between the fluid and solid is the sum of $\left(b\right)=1, 2, 3, \dots , m$ surface components and the volumetric transport coefficients, ${\dot{m}}_{gs\left(b\right)}$.

Concerning the change in the porosity with time, ${\partial }_t\epsilon =?$, there are a few strategies to solve this problem [69]. Some researchers use quite an independent evolution equation for the porosity. In the present approach, since composite materials can be strongly deformed, the change in the porosity can be described by an additional “closure” equation, i.e., Equation (44). Additionally, it was assumed that the pore surface density, $a$, and the tortuosity, $\tau$, are constant.

4.1.2. Evolution of Species in Reacting Gases

The composition of the homogeneous and heterogeneous evolutions of $N_R+m$ species within the reacting mixture is given in Equations (1) and (18). However, these equations cannot be simply added, since Equation (18) must be multiplied by the surface density, $a$. Therefore, the set of governing equations must also contain the homogeneous and heterogeneous (catalytic) evolution for $N_s-1$ species, where the rate of the increase/decrease in the species content mainly depends on chemical reactions:

$$\begin{gathered} {\partial }_t\left(\epsilon \rho Y_i+a{\epsilon }_s{\rho }_{cat}{Y}_{\left(i-m\right)}\right)+\mathrm{div}\left(\epsilon \rho Y_i{v}_g+a{\epsilon }_s{\rho }_{cat}{Y}_{\left(i-m\right)}s\right)= \\ =\mathrm{div}\left(J_i+a{j}_{\left(i-m\right)}\right)+\epsilon {\dot{\omega }}_i{W}_i+{\epsilon }_{a}a{\dot{\omega }}_{\left(i-m\right)}{W}_{\left(i-m\right)}+{\dot{m}}_{gs \left(i\right) }\left(Y_i-Y_{\left(i-m\right)}\right) \\ \mathrm{for} i=1,2,3,\dots ,N_s+m-1 \end{gathered}$$

(27)

where $Y_i$ (kg/kg) is the mass fraction of species $i=O_2,N_2,C{H}_4, H_2,\dots ,N_S$, which takes into account not only the bulk reaction but also the surface catalytic reaction. $J_i$ and $J_{\left(i-m\right)}$ are the Dixon–Levis diffusive flux bulk and surface transport, and $W_i$ is the molar density of each species. In addition, in Equation (27), the rate of the production of every species, ${\dot{\omega }}_i$, is related to the number of moles of the mixture.

4.2. Momentum Balance

A two-component continuum particle possesses two momentum balance equations, one to determine the gas velocity $v_g$ and the second to determine the solid velocity $v_s$:

$${\partial }_t\left(\epsilon {\rho }_g{v}_g\right)+\mathrm{div}\left(\epsilon {\rho }_g{v}_g⨂v_g\right)=\mathrm{div}\left({\tau }_g-p_{g}I+{\sigma }_g\right)+b_{gs}\left(v_g-v_s\right)$$

(28)

$${\partial }_t\left(\left(1-\epsilon \right){\rho }_s{v}_s\right)+\mathrm{div}\left(\left(1-\epsilon \right){\rho }_s{v}_s⨂v_s\right)=-\mathrm{div}\left({\sigma }_s-p_{g}I\right)-b_{gs}\left(v_g-v_s\right)$$

(29)

In the above (six-scalar) Equations (28) and (29), $\epsilon$ is the porosity, ${\rho }_g and {\rho }_s$ are the densities of the mixture of gases and solid composites, $v_g and v_s$ are the velocity of the gas mixture and the solid, ${\tau }_g$ is the total viscous stress tensor in the gas mixture, $p_g$ is the pressure of the gas mixture calculated based on the Dalton assumption about the summation of partial pressures of gas components, $b_{gs}\left(v_g-v_s\right)$ is the fluid–solid momentum interaction determined by the local coefficient of interaction $b_{gs}$, ${\sigma }_s$ is the partial stress tensor at the solid composite, and, finally, ${\sigma }_{terz,s}={\sigma }_s-p_{g}I$ is the so-called effective Terzaghi stress tensor. The direct presence of the fluid pressure, $p_g$, in the solid momentum balance is due to some freedom of fluid penetration inside the solid; however. such a directly symmetrical influence of the solid on the momentum flux tensor does not yet exist. Therefore, the introduction of the term ${\sigma }_g$, which is a recoverable, elastic-like stress tensor in the gas component, comes from the fluid–solid interaction. By analogy, the term $t_{terz,g}={\tau }_g-p_{g}I+{\sigma }_g$ is called the effective Terzaghi stress tensor in the bulk gas. Note that in Equations (28) and (29), classical body forces were omitted.

There is a special momentum balances case when a solid is deformable but stationary ($v_s=0$). leading to a special form of Equations (28) and (29). The momentum balance in the gas mixture of the continuum is given as follows:

$${\partial }_t\left(\epsilon {\rho }_g{v}_g\right)+\mathrm{div}\left(\epsilon {\rho }_g{v}_g⨂v_g\right)=\mathrm{div}\left({\tau }_g-p_{g}I\right)+f_{D-F}$$

(30)

where the resultant force $f_{D-F}$ is called the Darcy–Forchheimer slip resistance force (sometimes called “the pressure drop”) (Moghaddam and Jamiolahmady [70]). It is given as:

$$f_{D-F}=\frac{\mu }{K_1} \epsilon v_g+\frac{\mu }{K_2} {\epsilon }^2\left|v_g\right|v_g$$

(31)

where $\mu$ is the bulk shear viscosity of the mixture, and $K_1 and K_2$ are the Darcy–Forchheimer slip resistance coefficients (Moghaddam and Jamiolahmady [70], Sobieski and Trykozko [17]).

In the literature, together with fluid momentum balance, there is also a fluid porosity balance, which is written as an evolution equation for porosity $\epsilon$:

$${\partial }_t\left({\rho }_g\epsilon \right)+\mathrm{div}\left({\rho }_g{j}_D\right)=J$$

(32)

where ${\rho }_g$ is the fluid density, $\epsilon$ is the porosity, $j_D$ (a vector quantity) is the Darcy flux, or more specifically, the flux of the pore fluid. The first term in Equation (31) describes the changes in the fluid stored in the control volume, and the second term describes the net fluid flux across the control volume faces. In the above, $J$ is a fluid source or sink, which may be either at discrete locations or distributed in an arbitrary fashion throughout the domain. It is an open question how to derive the Darcy flux, $j_D$, because some force equilibrium conditions apply to the fluid as well as the porous medium. Constitutive closure is manifested through a relationship between the fluid flux and the forces driving the flux. For instance, for nearly all geological applications, the appropriate constitutive relationship is Darcy law, which can be expressed as (Krakowska et al. [18]):

$$j_D=-k\left(\frac{1}{{\mu }_g}\right)\left[\mathrm{grad} p_g+{\varrho }_{g}g \mathrm{grad} \psi +p_g \mathrm{grad}\left(I_{\epsilon }\right) \right]$$

(33)

where $k=k_{ml}{e}_m⨂e_l$ is a second-order permeability tensor, ${\mu }_g$ is the fluid dynamic viscosity, $g$ is the gravitational acceleration, $\psi$ is the gravitational potential (elevation above an arbitrary datum), and $I_{\epsilon }=\mathrm{tr} \left(\epsilon \right)$ is the trace of the solid deformation tensor. Darcy’s law shows that the porosity change is a function of the flow fluxes driven by gradients in terms of pressure energy, deformation energy, and elevation energy. In the undertaken approach, it was assumed that the porosity only changes due to deformation.

The static momentum balance in a porous solid (Equation (29)) reduces to the form:

$$\mathrm{div}\left({\sigma }_s-p_{g}I\right)-f_{D-F}=0$$

(34)

4.3. Three-Field Energy Balance

According to the concept of three temperature fields, every particle control volume needs three independent energy balance equations for all the unknowns:

4.3.1. Energy Balance in the Gas Constituent

The energy balance in the gas constituent is related to the internal energy balance of the gas mixture, $ϵ=c_p{T}_g$, in the following form:

$${\partial }_t\left(\epsilon {\rho }_gϵ\right)+\mathrm{div}\left(\epsilon \rho ϵv_g\right)=\mathrm{div} \left(\epsilon q_g\right)+h_{cat}\left(T_{cat}-T_g\right)-\epsilon \left({{\displaystyle \sum }}_{k=1}^{N_R}{\dot{\omega }}_k{h}_k\right)$$

(35)

In the three-temperature model, $T_{cat}$, $T_s$, and $T_g$, there are three interpenetrating temperatures in the continuum volume that cannot be simply mixed and equilibrated. Therefore, a “volumetric” heat exchange coefficient, $h_{cat}$, needs to be introduced. Note that in the first step, the contribution due to the surface catalytic temperature can be represented as a surface source: $h_{cat}\left(T_{cat}-T_s\right)\cong {\displaystyle \sum }_{\left(b\right)=1}^m{\dot{\omega }}_{\left(b\right)}{W}_{\left(b\right)}{h}_{\left(b\right)}$.

Equation (35) also contains classical parts that are additionally equilibrated with an internal energy exchange between the fluid and solid. There are two types of sources, one due to chemical reactions in the bulk gas, $N_R$, and a second due to the catalytic reactions on the surface, $N_c$. In this equation, the mixture density and heat capacity can be defined as:

$${\rho }_g={{\displaystyle \sum }}_{i=1}^{N_S}{Y}_i{\rho }_i, c_p={{\displaystyle \sum }}_{i=1}^{N_S}{Y}_i{c}_{p,i}$$

(36)

The rate of the bulk gas reaction is denoted as ${\dot{\omega }}_k$, and the enthalpy of the reaction is given by $h_k$. Similarly, the rate of the catalytic reaction at the surface is denoted as ${\dot{\omega }}_{\left(b\right)}$, and the enthalpy of the catalytic reaction is given by $h_{\left(b\right)}$.

4.3.2. Energy Balance in the Solid Constituent

The energy balance in the solid constituent takes into account heat diffusion flow, volumetric heat exchange, and energy radiation to the environment. It is given as:

$${\partial }_t\left(\left(1-\epsilon \right){\rho }_s{ϵ}_s\right)=\mathrm{div} \left(\left(1-\epsilon \right)q_s\right)-h_{cat}\left(T_{cat}-T_s\right)-{\dot{Q}}_{rad}$$

(37)

where the solid energy density is ${ϵ}_s=c_p{T}_s$, $q_s$ is the heat flux in the solid, and the main source of energy is emitted through radiation, ${\dot{Q}}_{rad}$ [J s⁻¹ m⁻²].

4.3.3. Energy Balance in the Catalytic Layer in Terms of Temperature

In the multi-field concept of a continuum, the distribution of the catalytic temperature, $T_{cat}$, cannot be established directly from Equation (24) due to the too many contact faces between the gas and the surface of the porous media. As a result, some averaging technique needs to be applied. Assuming that the surface density is constant, $a=const$, Equation (24) can be multiplied by the surface density to level-off the dimensions and spear the field of the catalytic temperature. Under these assumptions, it transforms into the following surface-like energy equation:

$$\begin{gathered} {\partial }_t\left(a{\rho }_{cat}{c}_{cat}{T}_{cat}\right)+\mathrm{div}\left(a{\rho }_{cat}{c}_{cat}s\right)+\mathrm{div}\left(a{j}_{cat}\right)+h_{cat}\left(T_g-T_s\right)+ \\ +{\sigma }_{em}{\epsilon }_{em}\left(T_{cat}^4-T_{ref,c}^4\right)+{\displaystyle \sum }_{\left(b\right)=1}^m{\dot{\omega }}_{\left(b\right)}{W}_{\left(b\right)}{h}_{\left(b\right)}=0 \end{gathered}$$

(38)

where $s$ is the Stefan catalytic slip velocity, $c_{cat}$ denotes the specific heat capacity of the catalyst, and ${\rho }_{cat}$ is the density of the catalyst layer. The surface-like heat flux, $j_{cat}$, is given by an analogy to the Fourier conductivity law, i.e., $j_{cat}={\lambda }_{cat}\mathrm{grad}(T_{cat})$, and $h_{\left(b\right)}$ is the specific enthalpy of the catalytic species $b$, either in the gas phase or at the surface. In the radiation term, ${\sigma }_{em}$ is the Stefan–Boltzmann constant, ${\epsilon }_{em}$ is the temperature-dependent surface emissivity, and $T_{ref,c}$ is the reference temperature at which the surface radiates.

The volumetric heat exchange coefficient, $h_{cat}$, is the same as in Equations (35) and (37), allowing a direct connection of the three temperature fields.

5. Kinematic Relationships

Kinematic relationships are based on a continuum divided into particles containing the following basic unknowns: $N_s$ fields of the mass fraction of species $Y_i$, $N_s$ fields of the mass fraction of the surface species $Y_{\left(b\right)}$, the actual solid density, two fields of velocity $v_g$ and $v_s$, and three fields of temperature, $T_{cat}$, $T_s$, and $T_g$. Therefore, the kinematic relationships are as simple as possible. The most important parameter is the tensor of the rate of fluid deformation, which is defined as a symmetrical part of the velocity gradient:

$$d=\frac{1}{2} \left(\mathrm{grad} v_g+{\mathrm{grad}}^{\mathrm{T}} v_g\right)=d_{jk}{e}_j⨂e_k$$

(39)

where the first invariant is equal to $I_d=\mathrm{tr}\left(d\right)=\mathrm{div} v_g$. For solids, this kinematical relationship looks similar. The rate of solid deformation is given as:

$$\dot{\epsilon }=\frac{1}{2} \left(\mathrm{grad} v_s+{\mathrm{grad}}^{\mathrm{T}} v_s\right)={\dot{\epsilon }}_{jk}{e}_j⨂e_k$$

(40)

where the small deformation tensor, $\epsilon$, is based on solid displacement. The gradient of temperature fields is used in the fluid and solid sub-domains, for instance:

$$g_g=\mathrm{grad} T_g \mathrm{and} g_s=\mathrm{grad} T_s$$

(41)

where $g_g$ and $g_s$ are the temperature gradients in the fluid and solid, respectively.

6. Constitutive Relationship for the Interaction Coefficients

The presented model is based on three main coefficients, namely ${\dot{m}}_{gs\left(i\right)}$, $b_{gs}$, and $h_{cat}$ (appearing in Section 4.3), which are the volumetric mass transport, the volumetric momentum transport, and the volumetric heat exchange coefficient, respectively. All of them have not yet been determined in this paper.

6.1. Volumetric Mass Transport

The volumetric mass transport coefficient is responsible for the mass jump related with the Lewis jump of the concentration. It is correlated with the surface jump coefficient $L_{\left(i\right)}$ (Groppi et al. [11]), which is individually measured for every species’ mass concentration (fraction):

$${\dot{m}}_{gs\left(i\right)}=a{\phi }_{\left(i\right)}{L}_{\left(i\right)} \mathrm{for} i=\left(b\right)$$

(42)

where $a$ is the surface volumetric density, and ${\phi }_{\left(i\right)}$ is the site fraction (Equation (20)).

6.2. Volumetric Momentum Transport

The coefficient $b_{gs}$ links two momentum balance equations, namely Equations (28) and (29). It forms the interchange force, $b_{gs}\left(v_g-v_s\right)$, between the solid and fluid and vice versa. In this particular case, when the solid velocity is equal to zero ($v_s=0$), the expression $b_{gs}\left(v_g-v_s\right)$ is reduced to the Darcy–Forchheimer slip resistance force $f_{D-F}$. From the above considerations, $b_{gs}$ must be dependent on the following parameters:

$$b_{gs}=b_{gs}\left(\epsilon ,a,\tau ,\mu ,{\mu }^{″}, \dots ,\left|v_g\right|,p_g\right)$$

(43)

which includes the porosity, surface density, tortuosity, bulk and surface viscosity, slip speed $\left|v_g\right|$, and other surface contributions. The bulk viscosity $\left(\mu \right)$ and surface viscosity $({\mu }^{″})$ appear together as the Navier slip length $l_s=\mu /{\mu }^{″}$. Alternatively, the surface viscosity together with the surface density, $a$, form a physically measurable coefficient, $K_1=f\left({\mu }^{″},a\right)$, which is called the permeability coefficient in the literature (sometimes denoted by the letter $\kappa$). Physically, the permeability coefficient describes the interactions between the surface viscosity and the whole control volume. Historically, the first scientist who measured this coefficient was Henry Darcy (1856) [71].

Another factor is the tortuosity, $\tau$, which is strongly related to the slip velocity, $\left|v_g\right|$, via the second permeability coefficient, $K_2=f(\tau ,$ $\left|v_g\right|, a)$, also known as the inertial permeability or the Forchheimer coefficient. The importance of this parameter depends on the type of flow. In the case of slow flow in the middle of a a sandstone reservoir, the Forchheimer equation is usually not needed, but in the case of gas inflow into a gas production well, the velocity may be high enough to justify the use of the Forchheimer correction. Some carbonate porous bodies have many internal fractures, so the Darcy equation for multiphase flow is generalized in order to govern both flows, i.e., in the fractures and in the matrix (i.e., traditional porous rock). The irregular surface of the fracture walls and the high flow rate in the fractures may also justify the use of the Forchheimer equation.

Another key parameter that is helpful in describing volumetric momentum transport is the Duhem slip force, $f_r^{Du}=p_g{\mu }^0{v}_g/\left|v_g\right|$, which is a part of the gas resistivity (${\mu }^0$ is the Duhem surface static-like viscosity [54]).

The Klinkenberg correction to the Darcy permeability, $K_1$, is defined as:

$$K_1^{eff}=K_1\left(1+\frac{b}{p_g}\right)$$

(44)

where the collective influence of the Duhem resistivity, ${\mu }^0$, is now represented by the Klinkenberg parameter. Finally, the volumetric momentum transport is represented via the following formula:

$$b_{gs}=\frac{\mu }{K_1\left(1+b/p_b\right)} \epsilon +\frac{\mu }{K_2} {\epsilon }^2\left|v_g\right|$$

(45)

$K_1$ is measured in m². An additional important factor is the hydraulic conductivity coefficient, $K=K_1\rho g/\mu$, which is measured in m/s (see Sobieski and Trykozko [17]).

6.3. Volumetric Heat Transport

The key parameter governing the volumetric heat exchange is the coefficient $h_{cat}$, which was already used in Equations (35) and (37). This coefficient can be defined on the basis of the surface heat transport equation, Equation (24):

$$h_{cat}=a_{v}h \mathrm{and} a_v=169.4 PPC$$

(46)

where $h$ is the surface heat transfer coefficient (between the gas mixture and the solid), $a_v$ is the cross-sectional area related to the unit volume of the porous material, and $PPC$ (pores per centimeter) is the number of pores per centimeter in the length of the porous medium. Please note the difference between the convective heat transfer coefficient $h$ and the volumetric heat transfer coefficient $h_{cat}$(see Nomenclature).

7. Constitutive Relationship for Bulk Properties

For fluids, the equation governing viscous Newtonian gas mixtures is given by (Badur and Ochrymiuk [8]):

$${\tau }_g=-\frac{2}{3}\mu I_{d}I+2\mu d$$

(47)

where $\mu$ is the gas mixture viscosity. The Fourier-like diffusional flux of thermal energy is given by:

$$q_g={\lambda }_g{g}_g$$

(48)

The heat diffusion in the gas mixture is governed by the total diffusion coefficient, ${\lambda }_g$ (Badur and Ochrymiuk [8]), and the gradient of the gas temperature is given by $g_g=\mathrm{grad} T_g$.

For solids, the heat flow is enhanced, since a collective mode of heat transport appears in the Fourier relation:

$$q_s={\lambda }_{s,eff}{g}_s$$

(49)

The heat diffusion in the solid is enhanced via additional contributions due to the gas flow velocity $v_g$ and is given as (Sobhan and Peterson [72]):

$${\lambda }_{s,eff}={\lambda }_s+\left|v_g\right| d c_{p,g}{K}_1^{-1}$$

(50)

However, the radiation energy source is given by the classical Stefan–Boltzmann formula (Weisz [14]):

$${\dot{Q}}_{rad}=\left(1-\epsilon \right) \alpha k_B\left(T_s^4-T_0^4\right)\exp \left(-\frac{\psi x}{l}\right)$$

(51)

where $k_B$ is the Stefan–Boltzmann constant, $\alpha$ is the emissivity of the solid, $\psi =-\ln \epsilon$ is the extinction function, and $x/l$ is the dimensionless position of the radiating place. The temperature of the environment was determined to be $T_0={288}^o \mathrm{K}$.

8. Terzaghi Concept of Effective Stresses in Porous Solids

The concept of mechanical coupling between fluids and porous solids has been recognized for over a century, as the first significant concept of FSIs combining models appeared in the 1920s. In Europe, Karl Terzaghi developed two crucial ideas, namely the notions of effective stress and the diffusion of fluid pressure due to flow. The classical Terzaghi concept of effective stress can be stated as (Kubik and Cieszko [71]):

$$\sigma ={\sigma }^{\prime }+\alpha p_{g}I=\left({\sigma }_{jk}^{\prime }+\alpha p_g{\delta }_{jk}\right)e_j⨂e_k$$

(52)

which is a constitutive equation of porous material. It proposes to decompose the elastic formulae of the stress tensor to account for the joint properties of the solid and surrounding fluid. Then, if:

$${\sigma }_{jk}^{\prime }=2G{\epsilon }_{jk}+K_s{\epsilon }_{mm}{\delta }_{jk}$$

(53)

is the classical Hooke law for a linearly elastic material, a change in the volumetric effective stress is given by three parameters: $G, K_s, and p_g$. In Equation (52), the correcting coefficient, $\alpha ,$ takes into account the influence of “pore rigidity” in the Biot form (Kubik and Cieszko [71]) $\alpha =1-K/K_s$ (where $K$ is the bulk modulus of the porous medium and $K_s$ is the bulk modulus of the solid grains). The change in the fluid pressure is denoted by $p$. The equation accurately describes the behavior of rocks under laboratory conditions. As the porosity approaches zero, $K$ gets closer to $K_s$, so the influence of the pore pressure on the effective stress vanishes as expected [17].

Under the assumption that the deformation of solid continua is “additive”, the linearized strain tensor $\epsilon =\frac{1}{2} \left(\mathrm{grad} u+{\mathrm{grad}}^{\mathrm{T}}u\right)={\epsilon }_{kl}{e}_k⨂e_l$ can be composed as a sum of particular independent contributions, such as elastic, porous, thermal, and gas sorption:

$$\begin{gathered} \epsilon ={\epsilon }^{el}+{\epsilon }^{po}+{\epsilon }^{th}+{\epsilon }^{so}= \\ ={\epsilon }^{el}+\frac{\alpha }{K} \left(p-p_0\right)I+\beta \left(T-T_0\right)I+{\displaystyle \sum }_{b=1}^m{\beta }_{\left(b\right)}\left(Y_{\left(b\right)}-Y_{0\left(b\right)}\right)I \end{gathered}$$

(54)

This means that elastic, porous, thermal, or gas sorption stresses cannot be indicated separately. Therefore, it was assumed that there are the following constitutive relations for isotropic materials between the total stress and total deformation tensors:

$$\begin{gathered} \sigma =2G\epsilon +2G\frac{\nu }{1-2\nu }tr\left(\epsilon \right)I+\alpha \left(p-p_0\right)I+ \\ +2G\frac{1-\nu }{1-2\nu }\left[{\beta }_T\left(T-T_0\right)+{\displaystyle \sum }_{\left(b\right)=1}^m{\beta }_{\left(b\right)}(Y_{\left(b\right)}-Y_{0\left(b\right)}\right]I \end{gathered}$$

(55)

This is in a form that is convenient for direct use in Equation (6) as $\mathrm{div} \sigma$. In the above, $G and v$ are the Kirchhoff and Poisson coefficients for solid (non-porous) materials, $\alpha$ is the Biot coefficient $\alpha =1-K/K_s$ (Kubik and Cieszko [71]) (where $K$ is the bulk modulus of the porous medium and $K_s$ is the bulk modulus of the solid grains), ${\beta }_T$ is the thermal expansion coefficient for porous material, and ${\beta }_{\left(b\right)}$ is the coefficient of gas sorption. The reference pressure, temperature, and mass fraction, i.e., $p_0, T_0, and Y_{0\left(b\right)}$, should be known beforehand.

9. Terzaghi Concept of Effective Stresses in Fluids

The physical closure for ${\sigma }_g$ needs to be discussed. This additional gas transport flux is related to the simplest phenomena of solid-body-like deformations of pores:

$${\sigma }_g=a^{-2}{\rho }_g{\left(s\right)}^2\mathrm{grad} I_{\epsilon }⨂\mathrm{grad} I_{\epsilon }$$

(56)

The conclusion from this is that the effective stresses in a fluid are not in the form of a spherical tensor (which is unusual in fluid constitutions), and it consists of shearing components such as ${\sigma }_{xy}$. In above equation, $a$ is the surface density, ${\rho }_g$ is the gas moisture density, $s$ is the Stefan slip velocity, and $I_{\epsilon }=\mathrm{div} u$ is the first invariant of the solid deformation tensor.

10. Terzaghi Effective Stresses and a Concept of Material Effort

Elasto-plastic and visco-elastic strains occur in the composite material of a burner. Different modes of plastic yielding are observed in composites during constant strain tests. Whether a composite experiences plastic yielding or suffers brittle failure depends on the stress state, temperature, and chemical environment. Plastic deformation can occur as combination of a shear and a volumetric strain, as seen in the study by Kubik and Cieszko [71] and Shi et al. [40]. Visco-elastic behavior (also called creep or time-dependent deformation in composite mechanics) can be very important in in-service operation during start-up and shutdown. Then, the constitutive relationship in Equation (52) must be extended to include the additional contribution of plastic and viscous strains.

A key challenge with elasto-plasticity (and elasto-viscosity) is the specification of a material effort by defining a yield surface that describes the boundary between the elastic and plastic state of a material. This problem has been well approached in terms of the principal stresses, ${\sigma }_1,{\sigma }_2, and {\sigma }_3$, of the stress tensor, $\sigma ={\sigma }_1{t}_1⨂t_1+{\sigma }_2{t}_2⨂t_2+{\sigma }_3{t}_3⨂t_3$, which are the largest, intermediate, and lowest principal stresses and are related to the principal axes $t_1,t_2, and t_3$. Experimental results take into account the total strength, so in order to obtain the effective Terzaghi stresses, ${\sigma }^{\prime }$, carried by the skeleton, the pore pressure needs to be subtracted from the total stress tensor, $\sigma$. The basic approach is to treat the yield surface as a Mohr–Coulomb failure envelope, which predicts yielding when the effective principal stresses deviate sufficiently from a hydrostatic state. A more comprehensive approach is provided by the critical state theory, which postulates a yield surface dependency on the effective stresses, strain state, and load history.

Wound ceramic matrix composites (CMC) are built from ceramic fibers embedded in a ceramic porous matrix [40]. They have many advantages over classical materials due to their high oxidation resistance, graceful failure behavior, thermal shock resistance, and sufficient strength at elevated temperature. These properties exceed all other materials used in porous burner technology. The concept of material effort within such strongly anisotropic porous materials is very complex and is based on a general hypothesis proposed by von Mises (1928) [40]:

$$f={\sigma }_{ij}^{\prime }{F}_{ij}+{\sigma }_{ij}^{\prime }{\sigma }_{kl}^{\prime }{F}_{ijkl}-1=0 \mathrm{for} \left(i,j,k,l=x,y,z\right)$$

(57)

This concept of material effort is based on effective Terzaghi stresses; therefore, the strength limits, $F_{ij}$ (second-order tensor) and $F_{ijkl}$ (fourth-order tensor), should be experimentally realized within the proper pressure conditions [40]. Analyzing the ceramic composite WHIPOX has resulted in the determination of five independent limits, $F_1,F_2,F_{11},F_{22}, and F_{66}$, and has reduced (57) to the well-known Tsai-Wu hypothesis. Unfortunately, their experiments concerned only atmospheric pressure, where $p_g=p_{atm}$.

11. The Boundary Conditions

In the most general case, a porous combustion chamber is in the form of a long cylinder (Groppi et al. [11]) with circular inlet ($z=0)$ and outlet ($z=l$) sections. Based on a given inlet temperature $T_{g,in}$, inlet mixture pressure $p_{g,in}$, outlet mixture pressure $p_{g,out}$, and mixture content, the mass flow rate at the inlet can be calculated as $\dot{m}=\iint \rho v_g·n dA$. Additionally, the radiative inlet conditions for the solid ($q_{s,z=0}·n=q_{s,rad}$) should be correctly predicted. At the initialization step, the temperature of the catalytic layer can be set to be equal to the gas temperature, i.e., $T_{cat,in}=T_{g,in}$, and the catalytic density can be set based on the reference, out–off reaction density ${\rho }_{cat,in}={\rho }_{cat,ref}$.

12. The Rate of Disaggregated Energy Production

The Bejan concept of disaggregated energy production [34] (similar to the Feidt concept [35,36,44]) for the case of catalytic combustion, includes the explanation of key assumptions, such the whole production of local entropy. It depends on the generation of entropy due to heat, friction, diffusion, and chemical contributions:

$${\dot{E}}_{gen}={\dot{E}}_{gen,heat}+{\dot{E}}_{gen, fric}+{\dot{E}}_{gen, diff}+{\dot{E}}_{gen, chem}$$

(58)

In the above, ${\dot{E}}_{gen}$ (J s⁻¹) is the whole-domain (integral) energy rate disaggregation, which is produced over the whole computational domain. In a cylindrical coordinate system, ${\dot{E}}_{gen}$ for a cylindrical domain is given as:

$${\dot{E}}_{gen}=\iint 2\pi r {\dot{ϵ}}_{gen} drdx$$

(59)

where $r and x$ are the radial and axial coordinates, respectively. The following components of disaggregated energy production were adopted from the literature:

-

Generation due to the heat diffusional flow:

$${\dot{ϵ}}_{gen,heat}=T_g^{-2} {\lambda }_g[{({\partial }_x{T}_g)}^2+\left({\partial }_r{T}_g{)}^2\right]+H{\left(T_s-T_g\right)}^2$$

(60)

-

Generation due to viscous bulk and surface friction:

$${\dot{ϵ}}_{gen,fric}=T_g^{-1}\mu \mathsf{\Phi }+f_{D-F}·s \mathrm{and} \mathsf{\Phi }=2d·d$$

(61)

-

Production due to species mixing:

$${\dot{ϵ}}_{gen,diff}={{\displaystyle \sum }}_{i=1}^{N_S}\rho D_{mi} J_i·J_i+{{\displaystyle \sum }}_{\left(b\right)=1}^m{\rho }_{cat}{D}_{m\left(b\right)} j_{\left(b\right)}·j_{\left(b\right)}$$

(62)

-

Production due to catalytic chemical reactions:

$${\dot{ϵ}}_{gen,chem}={{\displaystyle \sum }}_{i=1}^{N_S}{Y}_i{h}_i{\dot{\omega }}_i+{\displaystyle \sum }_{\left(b\right)=1}^m{\dot{\omega }}_{\left(b\right)}{W}_{\left(b\right)}{h}_{\left(b\right)}$$

(63)

The above relations can be computed after the simulation in a post-processing stage.

13. An Example of Heat Transfer Analysis Using Thermal FSI

The elements of the above multi-field model have been verified within the scope of different publications (Badur et al. [7], Badur and Ochrymiuk [8], Lewandowski et al. [73], Ochrymiuk [27,27], Ziółkowski et al. [43], Ziółkowski and Badur [33]). As an example, its application can be demonstrated with a simple thermal FSI benchmark example of the Stanton experiment (Stanton [74]), which was prepared in order to prove the Reynolds mechanism of thermal energy transport (Reynolds [56], Stanton and Pannell [75]). This is a simple case, where heat flux flows across the cylindrical metal wall from the hot water stream ($T_1=$ 39.6 °C) towards a co-flowing cold water stream ($t_1=$18 °C), as shown in Figure 5.

The total length of pipes was 67 cm, and the internal diameters were 13.9 mm (cold) and 15.7 mm (hot) (Stanton [74]). The mass flux of hot water was constant during the experiment and was equal to 148 g/s, whereas the mass flux of cold water varied from 148 g/s to 27 g/s. Note that the cold water inlet temperature was kept constant.

The aim of the Stanton experiment was to demonstrate the role of turbulent heat transfer. For the same mass flux of 148 g/s of hot and cold water, the result was symmetrical due to the conservation of energy, i.e., an increase in the temperature of the cold water by 4.23 °C ($t_1=18$°C → $t_2=22.23$°C) induced at temperature drop in the hot water by the same amount ($T_1=39.6$°C → $T_2=35.37$°C). This result was consistent between the experiment and simulation (see

Table 1. Comparison of the experimental and numerical outlet temperature, $t_2$.

Cooling Water Mass Flow [g/s]	Velocity of Cooling Water [cm/s]	Stanton ${\mathit{t}}_2$ [°C]	CFD ${\mathit{t}}_2$ [°C]
148	98	22.23	21.72
104	69	22.45	22.17
88	58.1	22.60	22.47
66	43.6	22.76	23.25
43	28.6	23.05	24.71
27	18	23.25	23.59

). The cold water mass flux reduction from 148 g/s to 27 g/s resulted in a systematic increase in the outlet temperature, $t_2$, from 22.23 °C to 23.25 °C and, simultaneously, an increase in $T_2$ from 35.2 °C to 37.4 °C.

It is worth mentioning that quite similar results have been numerically obtained by implementing the Reynolds–Stanton analogy (Reynolds [76]). It was calculated that during the change in the mass flow rate of the cold water (average velocity change from 24.4 to 4.5 cm/s), the total exchange of energy due to heating increased from $\dot{H}=$ 336 W to $\dot{H}=$ 390 W and harmonized with the formulae proposed by Reynolds and empirically proved by Stanton (in an original notation):

$$\dot{H}=\left(A+B\rho c\right)\left(T_1-t_1\right)$$

(64)

with $A=343$ and $B=334$, as calculated by Stanton and repeated in thermal FSI from an overall integration of the normal component of the heat vector:

$$\dot{H}=\iint q·n dA$$

(65)

It has been shown by numerical simulation that the heat flux (Equation (64)) depends directly on the length of the heat exchanging pipe; therefore, the Reynolds–Stanton formulae should be replaced by:

$$\dot{H}=l\left(518+499\rho c\right)\left(T_1-t_1\right)$$

(66)

where $l$ is the length of the heat exchange pipe. It is also important that the heat flux, $q=q^l+q^t$, was very sensitive to the turbulent heat flux modelling method, which was in agreement with Karcz and Badur [30] and Rup and Wais [77].

14. Conclusions

Everyone agrees that mathematical modelling represents a powerful method that is able to support the design process at every stage of the life cycle of catalytic porous combustors used in gas turbines. Therefore, a mathematical model that is robust and not time-consuming needs to be developed. In the paper, we proposed a two-stage model, the first dealing with a relatively exact description of the catalytic layer. Thus, the equations of mass, momentum, and energy conservation within the boundary of this catalytic layer were proposed and discussed. The basic novelty of this model lies in the detailed modeling of surface geometry and surface kinematics.

In the second stage of modeling, we proposed an averaged “two/three fields equation” for porous catalytic combustors. The concept is based on the volumetric (not surface, so far) coefficients of mass, momentum, and energy exchange.

However, the number of complex physical and chemical phenomena affecting the performance of a catalytic combustor is so large that the identification and application of certain approximations is inevitable. Therefore, in the numerical example, one possibility for identifying the parameters of the model was shown. The thermodynamical framework of the model was also opened to further research, i.e., the Bejan approach can be adopted.

The literature contains numerous models, the most popular being direct numerical simulation (DNS), large eddy simulation (LES), the lattice Boltzmann model (LBM), density functional theory (DFT), molecular dynamics (MD), and the Monte Carlo (MC) method. In opposition to these, in the present paper, we developed a model based on the concept of fluid–solid interactions. The proposed model has some advantages. It is based on the number of zones taken into consideration (fluid, solid, and catalytic layer). Depending on this, the proposed model can be described as a two-mass (Equations (25) and (26)), two-momentum (Equations (28) and (29)), and three-temperature model (Equations (35), (37), and (38)). Finally, the key aspects of the different phenomena occurring at the catalyst section were assessed by comparing simulation results.