# A phase-splitting approach to describe macroscopically non-equilibrium transport in porous media

Michel Quintard, Yohan Davit, Brian D. Wood

► To cite this version:

Michel Quintard, Yohan Davit, Brian D. Wood. A phase-splitting approach to describe macroscopically non-equilibrium transport in porous media. 9th European Thermal Sciences Conference (Eurotherm 2024), Jun 2024, Lake Bled, Slovenia. pp.012003, 10.1088/1742-6596/2766/1/012003 hal-04605649

HAL Id: hal-04605649

https://ut3-toulouseinp.hal.science/hal-04605649

Submitted on 7 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A phase-splitting approach to describe macroscopically non-equilibrium transport in porous media

Michel Quintard, Yohan Davit, Brian D. Wood

► To cite this version:

Michel Quintard, Yohan Davit, Brian D. Wood. A phase-splitting approach to describe macroscopically non-equilibrium transport in porous media. 9th European Thermal Sciences Conference (Eurotherm 2024), Jun 2024, Lake Bled, Slovenia. pp.012003, 10.1088/1742-6596/2766/1/012003. hal-04605649

HAL Id: hal-04605649

https://ut3-toulouseinp.hal.science/hal-04605649v1

Submitted on 7 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# A phase-splitting approach to describe macroscopically non-equilibrium transport in porous media

To cite this article: M. Quintard et al 2024 J. Phys.: Conf. Ser. 2766 012003

View the [article online](#) for updates and enhancements.

You may also like

- Noise and thermal depinning of Wigner crystals

C Reichhardt and C J O Reichhardt

- Traffic flow velocity disturbance characteristics and control strategy at the bottleneck of expressway

Jun-Wei Zeng, , Yong-Sheng Qian et al.

- Effect of violating the traffic light rule in the Biham-Middleton-Levine traffic flow model

Zhong-Jun Ding, Rui Jiang, Ming Li et al.

<figure><img src="image_2.png" /></figure>

The Electrochemical Society
Advancing solid state & electrochemical science & technology

DISCOVER
how sustainability
intersects with
electrochemistry & solid
state science research

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>
PAPER • OPEN ACCESS

# A phase-splitting approach to describe macroscopically non-equilibrium transport in porous media

To cite this article: M. Quintard et al 2024 J. Phys.: Conf. Ser. 2766 012003

View the [article online](#) for updates and enhancements.

You may also like

- Noise and thermal depinning of Wigner crystals

C Reichhardt and C J O Reichhardt

- Traffic flow velocity disturbance characteristics and control strategy at the bottleneck of expressway

Jun-Wei Zeng, , Yong-Sheng Qian et al.

- Effect of violating the traffic light rule in the Biham-Middleton-Levine traffic flow model

Zhong-Jun Ding, Rui Jiang, Ming Li et al.

<figure><img src="image_2.png" /></figure>

The Electrochemical Society
Advancing solid state & electrochemical science & technology

DISCOVER
how sustainability
intersects with
electrochemistry & solid
state science research

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>
# A phase-splitting approach to describe macroscopically non-equilibrium transport in porous media

M. Quintard¹, Y. Davit¹ and B. Wood²

¹Institut de Mécanique des Fluides de Toulouse, allée Soula, 31400 Toulouse - France

²Chemical, Biological, and Environmental Engineering 2047 Kelley Engineering Center
Corvallis, OR 97331 - USA

E-mail: michel.quintard@toulouse-inp.fr

**Abstract.** Classical macro-scale dispersion equation cannot correctly represent non-local, non-equilibrium effects, and alternatives had to be proposed to overcome these discrepancies. Two-equation models have been widely used and proved useful to incorporate, to some extent, local non-equilibrium effects. More accurate descriptions may be obtained through the introduction of additional complex terms in the one-equation dispersion equation, e.g. spatial and time convolutions. An alternative more practical route rely on the introduction of a N-equations description. For instance, instead of a two-temperature model, transfers in the solid and/or fluid phase maybe represented by several equations, each one defined by some splitting framework. Splitting of the solid phase has been proposed in the past (multi-rate models). Splitting of the flowing phase has a rustic illustration in the case of the Coats and Smith model, i.e., "flowing" and "stagnant" fluid equations. In this paper, we explore further the possibilities offered by N-equations splitting of the flowing phase. The N-equations macro-scale model is developed applying an averaging technique to N-pseudo phases obtained by splitting the flowing phase based on different criteria, e.g. the velocity histogram. Most important model properties are the matrices of dispersion tensors and exchange coefficients, which are provided by a system of specific "closure" problems. To illustrate this approach, the methodology is applied to Taylor's dispersion problem, in which non-local effects are produced by pore-scale spatially distributed inlet conditions. The proposed N-equations description was compared to estimates of the N-temperatures obtained from direct numerical simulations. The N-equations model was able to reproduce accurately the dynamic of pore-scale computations, for such reputedly non-homogenizable problems, with a relatively small number of equations.

## 1. Introduction

Classical dispersion equation is widely used for describing heat transfer and mass transport in porous media ([1]). Its emergence as an accurate averaged model requires that several time-scale and length-scale constraints must be fulfilled. The classical example is the Taylor dispersion equation describing the transport of a scalar in a tube ([2]): a parabolic equation with a dispersion coefficient scaling as Pe² is obtained as an asymptotic behavior, i.e., when lateral diffusion has smoothed the radial spatial variations. Abnormal behavior, often referred to as macro-scale non-fickian behavior, has been observed in many applications, and this triggered the development of modified macro-scale models ([3, 4]). Various mechanisms make the homogenization of a parabolic equations very difficult, thus leading to more complex models

<figure><img src="image_5.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
# A phase-splitting approach to describe macroscopically non-equilibrium transport in porous media

M. Quintard¹, Y. Davit¹ and B. Wood²

¹Institut de Mécanique des Fluides de Toulouse, allée Soula, 31400 Toulouse - France

²Chemical, Biological, and Environmental Engineering 2047 Kelley Engineering Center
Corvallis, OR 97331 - USA

E-mail: michel.quintard@toulouse-inp.fr

**Abstract.** Classical macro-scale dispersion equation cannot correctly represent non-local, non-equilibrium effects, and alternatives had to be proposed to overcome these discrepancies. Two-equation models have been widely used and proved useful to incorporate, to some extent, local non-equilibrium effects. More accurate descriptions may be obtained through the introduction of additional complex terms in the one-equation dispersion equation, e.g. spatial and time convolutions. An alternative more practical route rely on the introduction of a N-equations description. For instance, instead of a two-temperature model, transfers in the solid and/or fluid phase maybe represented by several equations, each one defined by some splitting framework. Splitting of the solid phase has been proposed in the past (multi-rate models). Splitting of the flowing phase has a rustic illustration in the case of the Coats and Smith model, i.e., "flowing" and "stagnant" fluid equations. In this paper, we explore further the possibilities offered by N-equations splitting of the flowing phase. The N-equations macro-scale model is developed applying an averaging technique to N-pseudo phases obtained by splitting the flowing phase based on different criteria, e.g. the velocity histogram. Most important model properties are the matrices of dispersion tensors and exchange coefficients, which are provided by a system of specific "closure" problems. To illustrate this approach, the methodology is applied to Taylor's dispersion problem, in which non-local effects are produced by pore-scale spatially distributed inlet conditions. The proposed N-equations description was compared to estimates of the N-temperatures obtained from direct numerical simulations. The N-equations model was able to reproduce accurately the dynamic of pore-scale computations, for such reputedly non-homogenizable problems, with a relatively small number of equations.

## 1. Introduction

Classical dispersion equation is widely used for describing heat transfer and mass transport in porous media ([1]). Its emergence as an accurate averaged model requires that several time-scale and length-scale constraints must be fulfilled. The classical example is the Taylor dispersion equation describing the transport of a scalar in a tube ([2]): a parabolic equation with a dispersion coefficient scaling as Pe² is obtained as an asymptotic behavior, i.e., when lateral diffusion has smoothed the radial spatial variations. Abnormal behavior, often referred to as macro-scale non-fickian behavior, has been observed in many applications, and this triggered the development of modified macro-scale models ([3, 4]). Various mechanisms make the homogenization of a parabolic equations very difficult, thus leading to more complex models

<figure><img src="image_5.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
than the classical dispersion equation: spatial heterogeneities, vicinity of boundaries, spatially distributed initial conditions, preferential paths, etc. Several models were proposed to handle these difficulties [5]. Considering a two-phase system, its behavior may obviously be modeled through Direct Numerical Simulation (DNS), with the requirement of large computational resources. It is tempting to represent a macro-scale 1D transport problem through such a structure by a single equation model, i.e., a 1-equation model. The simplest 1-Eq models, like the Taylor dispersion equation, may fail to reproduce complex behaviors, as discussed above. Maintaining a 1-Eq description, while trying to reproduce as much as possible of the complex features requires to add complicated terms, e.g. spatially and time dependent effective properties, or adopt mathematical forms different from the classical parabolic equation. An alternative route is the introduction of N-equation models. This was particularly used, for instance, to describe local non-equilibrium heat transfer in two-phase systems, were the averaged temperature of one of the phase is lagging behind the other. More equations are needed, but non-equilibrium or abnormal effects are reproduced with a minimum number of additional effective properties, that do not require the complexity of 1-Eq models under such non-equilibrium conditions, as illustrated in [6].

In these examples, the two-phase structure motivates the division into different macro-scale equations. However, various spatial and time scales, the principal factor for abnormal behavior, may be encountered within a single phase itself. Indeed, the Coats and Smith model [7] is a typical and historical example of the use of 2-Eq model in such cases. In this model, the flowing phase was split into the mobile part, exchanging with the immobile part (e.g. dead-end pores). This was also suggested in the case of structured media used in distillation columns. Such media are made of corrugated sheets, with families of grooves having different directions, possibly triggering macro-scale different velocity orientations. Phase splitting according to this topological feature proved to offer an accurate description of such mechanisms [8, 9, 10]. This concept of phase-splitting was investigated from a general upscaling point of view in [11]. The approach was applied successfully to transport in nodular systems, with the splitting occurring in the solid phase, which was only subject to diffusion. In this paper, we explore further the concept by considering a splitting of the liquid phase, and abnormal effects produced by spatially distributed boundary or initial conditions. To illustrate the phase-splitting concept, we will consider a 2D flow between two planes distant of 2h. In the next section, we describe the pore-scale problem and the major upscaling steps leading to a N-Equation model. Then, we apply this framework to a 2D case, with spatially distributed entrance conditions. This choice emphasizes the ability of N-eqs models to capture the macro-scale behavior, without the need of complex expressions for the effective properties, or complicated partial derivatives terms.

## 2. Theoretical Framework

We consider the transport of a scalar quantity (e.g. temperature $T$) by a fluid flowing between two-planes distant from $2h$ (figure 1). The dimensionless pore-scale equation is written as

$$ \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} = \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \qquad (1) $$

where the dimensionless $x$-component of velocity, $u$, is given by

$$ \text{uniform flow: } u = Pe \quad ; \quad \text{Poiseuille flow: } u = \frac{3}{2}Pe(1 - y^2) \qquad (2) $$

The Péclet number involves $h$ as the reference length, all dimensionless variables are defined in the classical way. Problem is completed by giving initial and boundary conditions such that

$$ T(x, y, t = 0) = T_0(x, y) ; \quad \frac{\partial T}{\partial y} = 0 \quad \text{at } y = -1, +1 ; \quad T(x = 0, y, t) = g(y) \qquad (3) $$
than the classical dispersion equation: spatial heterogeneities, vicinity of boundaries, spatially distributed initial conditions, preferential paths, etc. Several models were proposed to handle these difficulties [5]. Considering a two-phase system, its behavior may obviously be modeled through Direct Numerical Simulation (DNS), with the requirement of large computational resources. It is tempting to represent a macro-scale 1D transport problem through such a structure by a single equation model, i.e., a 1-equation model. The simplest 1-Eq models, like the Taylor dispersion equation, may fail to reproduce complex behaviors, as discussed above. Maintaining a 1-Eq description, while trying to reproduce as much as possible of the complex features requires to add complicated terms, e.g. spatially and time dependent effective properties, or adopt mathematical forms different from the classical parabolic equation. An alternative route is the introduction of N-equation models. This was particularly used, for instance, to describe local non-equilibrium heat transfer in two-phase systems, were the averaged temperature of one of the phase is lagging behind the other. More equations are needed, but non-equilibrium or abnormal effects are reproduced with a minimum number of additional effective properties, that do not require the complexity of 1-Eq models under such non-equilibrium conditions, as illustrated in [6].

In these examples, the two-phase structure motivates the division into different macro-scale equations. However, various spatial and time scales, the principal factor for abnormal behavior, may be encountered within a single phase itself. Indeed, the Coats and Smith model [7] is a typical and historical example of the use of 2-Eq model in such cases. In this model, the flowing phase was split into the mobile part, exchanging with the immobile part (e.g. dead-end pores). This was also suggested in the case of structured media used in distillation columns. Such media are made of corrugated sheets, with families of grooves having different directions, possibly triggering macro-scale different velocity orientations. Phase splitting according to this topological feature proved to offer an accurate description of such mechanisms [8, 9, 10]. This concept of phase-splitting was investigated from a general upscaling point of view in [11]. The approach was applied successfully to transport in nodular systems, with the splitting occurring in the solid phase, which was only subject to diffusion. In this paper, we explore further the concept by considering a splitting of the liquid phase, and abnormal effects produced by spatially distributed boundary or initial conditions. To illustrate the phase-splitting concept, we will consider a 2D flow between two planes distant of 2h. In the next section, we describe the pore-scale problem and the major upscaling steps leading to a N-Equation model. Then, we apply this framework to a 2D case, with spatially distributed entrance conditions. This choice emphasizes the ability of N-eqs models to capture the macro-scale behavior, without the need of complex expressions for the effective properties, or complicated partial derivatives terms.

## 2. Theoretical Framework

We consider the transport of a scalar quantity (e.g. temperature $T$) by a fluid flowing between two-planes distant from $2h$ (figure 1). The dimensionless pore-scale equation is written as

$$ \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} = \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \qquad (1) $$

where the dimensionless $x$-component of velocity, $u$, is given by

$$ \text{uniform flow: } u = Pe \quad ; \quad \text{Poiseuille flow: } u = \frac{3}{2}Pe(1 - y^2) \qquad (2) $$

The Péclet number involves $h$ as the reference length, all dimensionless variables are defined in the classical way. Problem is completed by giving initial and boundary conditions such that

$$ T(x, y, t = 0) = T_0(x, y) ; \quad \frac{\partial T}{\partial y} = 0 \quad \text{at } y = -1, +1 ; \quad T(x = 0, y, t) = g(y) \qquad (3) $$
<figure><img src="image_6.png" /><figcaption>Figure 1. Computation domain, example of N-phase splitting.</figcaption></figure>

completed with outlet convective conditions.

The phase splitting procedure involves decomposing the domain into N subdomains, as illustrated figure 1. Each decomposed domain, i, is defined by the following interval

$$i : [y_{i-1}, y_i] \qquad (4)$$

with $y_0 = -1$. These intervals may be chosen arbitrarily, while they must realize a partition of the entire domain. Volume fractions are obtained as $\phi_i = \frac{y_i - y_{i-1}}{2}$. The i-phase intrinsic averaged velocity and intrinsic averaged temperature are given by

$$\phi_i U_i = \frac{1}{2} \int_{y^{(i-1)}}^{y^{(i)}} u(y) dy \quad ; \quad \mathbb{T}_i = \langle T_i \rangle^{(i)} = \frac{1}{y_i - y_{i-1}} \int_{y_{i-1}}^{y_i} T_i dy \qquad (5)$$

in which $T_i$ is the restriction of T to the i-domain. Temperature deviations, $\tilde{T}_i$, are defined as

$$T_i = \mathbb{T}_i + \tilde{T}_i \qquad (6)$$

With all these definitions, the i-phase averaged equation is given by

$$\frac{\partial \mathbb{T}_i}{\partial t} + \left\langle u \frac{\partial T_i}{\partial x} \right\rangle^{(i)} = \left\langle \frac{\partial^2 T_i}{\partial x^2} + \frac{\partial^2 T_i}{\partial y^2} \right\rangle^{(i)} \qquad (7)$$

Using this into equation 1, one can write a governing equation for the temperature deviations as

$$\frac{\partial \tilde{T}_i}{\partial t} + u \frac{\partial T_i}{\partial x} - \left\langle u \frac{\partial T_i}{\partial x} \right\rangle^{(i)} = \frac{\partial^2 T_i}{\partial x^2} + \frac{\partial^2 T_i}{\partial y^2} - \left\langle \frac{\partial^2 T_i}{\partial x^2} + \frac{\partial^2 T_i}{\partial y^2} \right\rangle^{(i)} \qquad (8)$$

We follow the route discussed in [11] to propose an approximate solution of the coupled macro-scale and pore-scale equations using the following representations of the temperature deviations with respect to the pore-scale deviation equations source terms

$$\tilde{T}_i = (a_{ij} - \delta_{ij}) \mathbb{T}_j + b_{ij} \text{Pe} \frac{\partial \mathbb{T}_j}{\partial x} \qquad (9)$$

where the mapping matrices, a and b, are solution of the following local problems

**Problem I:**

$$\frac{\partial^2 a_{ij}}{\partial y^2} = \left\langle \frac{\partial^2 a_{ij}}{\partial y^2} \right\rangle^{(i)} \qquad (10)$$

with the following boundary conditions

$$a_{ij} = a_{(i+1)j} \quad ; \quad \left. \frac{da_{ij}}{dy} \right|_{y_i} = \left. \frac{da_{(i+1)j}}{dy} \right|_{y_i} \qquad (11)$$
<figure><img src="image_6.png" /><figcaption>Figure 1. Computation domain, example of N-phase splitting.</figcaption></figure>

completed with outlet convective conditions.

The phase splitting procedure involves decomposing the domain into N subdomains, as illustrated figure 1. Each decomposed domain, i, is defined by the following interval

$$i : [y_{i-1}, y_i] \qquad (4)$$

with $y_0 = -1$. These intervals may be chosen arbitrarily, while they must realize a partition of the entire domain. Volume fractions are obtained as $\phi_i = \frac{y_i - y_{i-1}}{2}$. The i-phase intrinsic averaged velocity and intrinsic averaged temperature are given by

$$\phi_i U_i = \frac{1}{2} \int_{y^{(i-1)}}^{y^{(i)}} u(y) dy \quad ; \quad \mathbb{T}_i = \langle T_i \rangle^{(i)} = \frac{1}{y_i - y_{i-1}} \int_{y_{i-1}}^{y_i} T_i dy \qquad (5)$$

in which $T_i$ is the restriction of T to the i-domain. Temperature deviations, $\tilde{T}_i$, are defined as

$$T_i = \mathbb{T}_i + \tilde{T}_i \qquad (6)$$

With all these definitions, the i-phase averaged equation is given by

$$\frac{\partial \mathbb{T}_i}{\partial t} + \left\langle u \frac{\partial T_i}{\partial x} \right\rangle^{(i)} = \left\langle \frac{\partial^2 T_i}{\partial x^2} + \frac{\partial^2 T_i}{\partial y^2} \right\rangle^{(i)} \qquad (7)$$

Using this into equation 1, one can write a governing equation for the temperature deviations as

$$\frac{\partial \tilde{T}_i}{\partial t} + u \frac{\partial T_i}{\partial x} - \left\langle u \frac{\partial T_i}{\partial x} \right\rangle^{(i)} = \frac{\partial^2 T_i}{\partial x^2} + \frac{\partial^2 T_i}{\partial y^2} - \left\langle \frac{\partial^2 T_i}{\partial x^2} + \frac{\partial^2 T_i}{\partial y^2} \right\rangle^{(i)} \qquad (8)$$

We follow the route discussed in [11] to propose an approximate solution of the coupled macro-scale and pore-scale equations using the following representations of the temperature deviations with respect to the pore-scale deviation equations source terms

$$\tilde{T}_i = (a_{ij} - \delta_{ij}) \mathbb{T}_j + b_{ij} \text{Pe} \frac{\partial \mathbb{T}_j}{\partial x} \qquad (9)$$

where the mapping matrices, a and b, are solution of the following local problems

**Problem I:**

$$\frac{\partial^2 a_{ij}}{\partial y^2} = \left\langle \frac{\partial^2 a_{ij}}{\partial y^2} \right\rangle^{(i)} \qquad (10)$$

with the following boundary conditions

$$a_{ij} = a_{(i+1)j} \quad ; \quad \left. \frac{da_{ij}}{dy} \right|_{y_i} = \left. \frac{da_{(i+1)j}}{dy} \right|_{y_i} \qquad (11)$$
at the interface between i and i + 1 domains. We also have

$$ \frac{da_{1j}}{dy} = 0 \quad \text{at } y = -1; \quad \frac{da_{Nj}}{dy} = 0 \quad \text{at } y = +1; \quad \langle a_{ij} \rangle^{(i)} = \delta_{ij} \qquad (12) $$

### Problem II:

$$ \frac{d^2 b_{ij}}{dy^2} = \frac{u_i}{Pe} a_{ij} - \frac{U_i}{Pe} \delta_{ij} - \phi_i^{-1} C_{ij} \qquad (13) $$

where $C_{ij} = \langle \frac{\tilde{u}_i}{Pe} a_{ij} \rangle - \langle \frac{d^2 b_{ij}}{dy^2} \rangle$. The boundary conditions are given by

$$ b_{ij} = b_{(i+1)j} ; \quad \left. \frac{db_{ij}}{dy} \right|_{y_i} = \left. \frac{db_{(i+1)j}}{dy} \right|_{y_i} \qquad (14) $$

at the interface between i and i + 1 domains. We also have

$$ \frac{db_{1j}}{dy} = 0 \quad \text{at } y = -1; \quad \frac{db_{Nj}}{dy} = 0 \quad \text{at } y = +1; \quad \langle b_{ij} \rangle^{(i)} = 0 \qquad (15) $$

## 2.1. N-equations model

With the above approximation of the coupled pore-scale/macro-scale problems, one can finally write the N-equations model under the general form

$$ \left( \begin{array}{c} \frac{\partial \phi_1 \mathbb{T}_1}{\partial t} \\ \vdots \\ \frac{\partial \phi_N \mathbb{T}_N}{\partial t} \end{array} \right) + \left( \begin{array}{c} \phi_1 U_1 \frac{\partial \mathbb{T}_1}{\partial x} \\ \vdots \\ \phi_N U_N \frac{\partial \mathbb{T}_N}{\partial x} \end{array} \right) + \text{Pe} \left( \begin{array}{c} C_{11} \cdots C_{1N} \\ \vdots \cdots \vdots \\ C_{N1} \cdots C_{NN} \end{array} \right) \left( \begin{array}{c} \frac{\partial \mathbb{T}_1}{\partial x} \\ \vdots \\ \frac{\partial \mathbb{T}_N}{\partial x} \end{array} \right) = \left( \begin{array}{c} A_{11} \cdots A_{1N} \\ \vdots \cdots \vdots \\ A_{N1} \cdots A_{NN} \end{array} \right) \left( \begin{array}{c} \frac{\partial^2 \mathbb{T}_1}{\partial x^2} \\ \vdots \\ \frac{\partial^2 \mathbb{T}_N}{\partial x^2} \end{array} \right) - \left( \begin{array}{c} H_{11} \cdots H_{1N} \\ \vdots \cdots \vdots \\ H_{N1} \cdots H_{NN} \end{array} \right) \left( \begin{array}{c} \mathbb{T}_1 \\ \vdots \\ \mathbb{T}_N \end{array} \right) \qquad (16) $$

with $C_{ij}$ already defined, and

$$ A_{ij} = \langle a_{ij} \rangle - \text{Pe}^2 \left\langle \frac{u_i}{Pe} b_{ij} \right\rangle = \phi_i \delta_{ij} - \text{Pe}^2 \left\langle \frac{\tilde{u}_i}{Pe} b_{ij} \right\rangle ; \quad H_{ij} = - \left\langle \frac{\partial^2 a_{ij}}{\partial y^2} \right\rangle \qquad (17) $$

Besides the additional convective terms, the effective properties involve a matrix of Taylor dispersion coefficients and a matrix of exchange terms, equation 17. In the context of this 2D problem, solution of the "closure" problems I and II can be obtained analytically. A Python code was built, using symbolic calculus, to construct automatically the effective property matrices, given a chosen domain decomposition. This is an important feature of the proposed theoretical approach that effective properties obey relatively simple PDEs, so the number of domains can be adapted rapidly based on specific features of the boundary and initial conditions of the problem. Illustrative examples are provided in the next section.

## 3. Examples for a three domain decomposition

In this example, we adopt a domain decomposition made of equal size domains with $N = 3$. To illustrate the challenge that can be taken by the proposed approach, we solve several problems with the following initial and boundary conditions

$$ \left\{ \begin{array}{ll} \text{case I} & T_0 = 0; \quad g(y) = \mathcal{H}(y - y_1) - \mathcal{H}(y - y_2) \\ \text{case II} & T_0 = 0; \quad g(y) = \mathcal{H}(y - y_1) - \mathcal{H}(y - y_2) \end{array} \right. \qquad (18) $$
at the interface between i and i + 1 domains. We also have

$$ \frac{da_{1j}}{dy} = 0 \quad \text{at } y = -1; \quad \frac{da_{Nj}}{dy} = 0 \quad \text{at } y = +1; \quad \langle a_{ij} \rangle^{(i)} = \delta_{ij} \qquad (12) $$

### Problem II:

$$ \frac{d^2 b_{ij}}{dy^2} = \frac{u_i}{Pe} a_{ij} - \frac{U_i}{Pe} \delta_{ij} - \phi_i^{-1} C_{ij} \qquad (13) $$

where $C_{ij} = \langle \frac{\tilde{u}_i}{Pe} a_{ij} \rangle - \langle \frac{d^2 b_{ij}}{dy^2} \rangle$. The boundary conditions are given by

$$ b_{ij} = b_{(i+1)j} ; \quad \left. \frac{db_{ij}}{dy} \right|_{y_i} = \left. \frac{db_{(i+1)j}}{dy} \right|_{y_i} \qquad (14) $$

at the interface between i and i + 1 domains. We also have

$$ \frac{db_{1j}}{dy} = 0 \quad \text{at } y = -1; \quad \frac{db_{Nj}}{dy} = 0 \quad \text{at } y = +1; \quad \langle b_{ij} \rangle^{(i)} = 0 \qquad (15) $$

## 2.1. N-equations model

With the above approximation of the coupled pore-scale/macro-scale problems, one can finally write the N-equations model under the general form

$$ \left( \begin{array}{c} \frac{\partial \phi_1 \mathbb{T}_1}{\partial t} \\ \vdots \\ \frac{\partial \phi_N \mathbb{T}_N}{\partial t} \end{array} \right) + \left( \begin{array}{c} \phi_1 U_1 \frac{\partial \mathbb{T}_1}{\partial x} \\ \vdots \\ \phi_N U_N \frac{\partial \mathbb{T}_N}{\partial x} \end{array} \right) + \text{Pe} \left( \begin{array}{c} C_{11} \cdots C_{1N} \\ \vdots \cdots \vdots \\ C_{N1} \cdots C_{NN} \end{array} \right) \left( \begin{array}{c} \frac{\partial \mathbb{T}_1}{\partial x} \\ \vdots \\ \frac{\partial \mathbb{T}_N}{\partial x} \end{array} \right) = \left( \begin{array}{c} A_{11} \cdots A_{1N} \\ \vdots \cdots \vdots \\ A_{N1} \cdots A_{NN} \end{array} \right) \left( \begin{array}{c} \frac{\partial^2 \mathbb{T}_1}{\partial x^2} \\ \vdots \\ \frac{\partial^2 \mathbb{T}_N}{\partial x^2} \end{array} \right) - \left( \begin{array}{c} H_{11} \cdots H_{1N} \\ \vdots \cdots \vdots \\ H_{N1} \cdots H_{NN} \end{array} \right) \left( \begin{array}{c} \mathbb{T}_1 \\ \vdots \\ \mathbb{T}_N \end{array} \right) \qquad (16) $$

with $C_{ij}$ already defined, and

$$ A_{ij} = \langle a_{ij} \rangle - \text{Pe}^2 \left\langle \frac{u_i}{Pe} b_{ij} \right\rangle = \phi_i \delta_{ij} - \text{Pe}^2 \left\langle \frac{\tilde{u}_i}{Pe} b_{ij} \right\rangle ; \quad H_{ij} = - \left\langle \frac{\partial^2 a_{ij}}{\partial y^2} \right\rangle \qquad (17) $$

Besides the additional convective terms, the effective properties involve a matrix of Taylor dispersion coefficients and a matrix of exchange terms, equation 17. In the context of this 2D problem, solution of the "closure" problems I and II can be obtained analytically. A Python code was built, using symbolic calculus, to construct automatically the effective property matrices, given a chosen domain decomposition. This is an important feature of the proposed theoretical approach that effective properties obey relatively simple PDEs, so the number of domains can be adapted rapidly based on specific features of the boundary and initial conditions of the problem. Illustrative examples are provided in the next section.

## 3. Examples for a three domain decomposition

In this example, we adopt a domain decomposition made of equal size domains with $N = 3$. To illustrate the challenge that can be taken by the proposed approach, we solve several problems with the following initial and boundary conditions

$$ \left\{ \begin{array}{ll} \text{case I} & T_0 = 0; \quad g(y) = \mathcal{H}(y - y_1) - \mathcal{H}(y - y_2) \\ \text{case II} & T_0 = 0; \quad g(y) = \mathcal{H}(y - y_1) - \mathcal{H}(y - y_2) \end{array} \right. \qquad (18) $$
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Comparison between DNS averaged temperatures (denoted DNS[i-1] for phase i) and macro-scale predictions (denoted Macro[i-1] for phase i) from the 3-Eq. model (Pe = 10, t = 0.5)</figcaption></figure>

where *H* is the Heaviside "function". The analysis will follow the following steps: (i) DNS: the pore-scale equations are solved numerically with the help of a OpenFoam® code, with usual numerical parameters for this very classical parabolic PDE, (ii) Computation of averaged temperature fields from DNS, (iii) Computation of effective properties, (iv) Numerical solution of the macro-scale equations, (v) Comparison between macro-scale and averaged temperature fields. For these tests, we have discarded the C-matrix in order to simplify the macro-scale solver (in this case, there is no coupling between the advection terms, and an operator splitting technique may be used for the numerical implementation). However, assessing the importance of such terms requires certainly further analysis.

### 3.1. *Results for case I, uniform flow, symmetric BC*

The inlet boundary condition features a characteristic length-scale on the same order as *h*, and this result in a complicated mixing process. This complex behavior is illustrated in figure 2, which represents the averaged and macro-scale temperature fields as a function of *x*, for a given Pe number at a given time. One sees that the macro-scale model predicts well the exchange of heat between the 3 layers near inlet, and also the long range mixing process.

### 3.2. *Results for case II, Poiseuille flow, asymmetric BC*

In this case, we test further the 3-Eq approach by imposing an asymmetric inlet boundary condition. One sees a rather complicated mixing process, especially in the region around *x* = 6. This is confirmed when looking at the DNS averaged fields as shown in figure 3. While the predictions of the macro-scale model (with only this limited number of domains) are less accurate than in the previous case, one sees a relatively good agreement, and a correct prediction of the intricate evolution of the three fields around *x* = 6.

## 4. Conclusions

In this paper, we investigated the ability of a N-Eq macro-scale description to represent correctly the impact of spatially dependent boundary conditions. The chosen BC characteristic length was such that classical Taylor dispersion 1-equation fails to capture the induced complex behavior. Attempts to overcome this discrepancy through a one equation model would require to *propagate*
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Comparison between DNS averaged temperatures (denoted DNS[i-1] for phase i) and macro-scale predictions (denoted Macro[i-1] for phase i) from the 3-Eq. model (Pe = 10, t = 0.5)</figcaption></figure>

where *H* is the Heaviside "function". The analysis will follow the following steps: (i) DNS: the pore-scale equations are solved numerically with the help of a OpenFoam® code, with usual numerical parameters for this very classical parabolic PDE, (ii) Computation of averaged temperature fields from DNS, (iii) Computation of effective properties, (iv) Numerical solution of the macro-scale equations, (v) Comparison between macro-scale and averaged temperature fields. For these tests, we have discarded the C-matrix in order to simplify the macro-scale solver (in this case, there is no coupling between the advection terms, and an operator splitting technique may be used for the numerical implementation). However, assessing the importance of such terms requires certainly further analysis.

### 3.1. *Results for case I, uniform flow, symmetric BC*

The inlet boundary condition features a characteristic length-scale on the same order as *h*, and this result in a complicated mixing process. This complex behavior is illustrated in figure 2, which represents the averaged and macro-scale temperature fields as a function of *x*, for a given Pe number at a given time. One sees that the macro-scale model predicts well the exchange of heat between the 3 layers near inlet, and also the long range mixing process.

### 3.2. *Results for case II, Poiseuille flow, asymmetric BC*

In this case, we test further the 3-Eq approach by imposing an asymmetric inlet boundary condition. One sees a rather complicated mixing process, especially in the region around *x* = 6. This is confirmed when looking at the DNS averaged fields as shown in figure 3. While the predictions of the macro-scale model (with only this limited number of domains) are less accurate than in the previous case, one sees a relatively good agreement, and a correct prediction of the intricate evolution of the three fields around *x* = 6.

## 4. Conclusions

In this paper, we investigated the ability of a N-Eq macro-scale description to represent correctly the impact of spatially dependent boundary conditions. The chosen BC characteristic length was such that classical Taylor dispersion 1-equation fails to capture the induced complex behavior. Attempts to overcome this discrepancy through a one equation model would require to *propagate*
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Comparison between DNS averaged temperatures (denoted DNS[i-1] for phase *i*) and macro-scale predictions (denoted Macro[i-1] for phase *i*) from the 3-Eq. model (Pe = 10, *t* = 0.7)</figcaption></figure>

the macro-scale PDE with various terms and complex (possibly spatially and time dependent) effective properties. Most of the time, additional terms would be specific to the pore-scale initial boundary value problem under investigation.

The proposed alternative of a N-eq description offers more flexibility. First, we can adapt the N-eq description in terms of *N*, and also the domain-decomposition topology, to characteristic length-scales of the BCs spatial distribution. Like one-equation descriptions, N-equation models can have complex mathematical structures (e.g. spatial and time dependent effective properties in convolution products). However, the choice of *N* and domain topology is already a powerful tool to approach the actual macro-scale behavior. Therefore, we tested here a first-order closure for the N-eq model which is *independent* of the actual IBVP under investigation, and gives simple expressions for the effective properties. Indeed, in the particular case used as an illustration in this paper, effective properties are matrices with constant coefficients obtained from analytical expressions (the use of symbolic programming allows to construct automatically these properties for any *N*, domain topology, and velocity field). Various tests showed that the macro-scale predictions are rather robust estimates of the actual behavior.

## References

[1] Bear J 1972 *Dynamics of Fluids in Porous Media* (New York: Elsevier)

[2] Taylor G I 1953 *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences* **219** 186–203

[3] Koch D L and Brady J F 1987 *J. Fluid Mech.* **180** 387–403

[4] Cushman H J 1997 *The Physics of Fluids in Hierarchical Porous Media: Angstroms to Miles* (The Netherlands: Academic Publishers)

[5] Davit Y and Quintard M 2015 *Handbook of Porous Media* ed Vafai K (Taylor & Francis) pp 245–320 3rd ed

[6] Davarzani H, Marcoux M and Quintard M 2010 *International Journal of Heat and Mass Transfer* **53** 1514–1528

[7] Coats K H and Smith B D 1964 *Society of Petroleum Engineers Journal March* 73–84

[8] Mahr B and Mewes D 2008 *AIChE Journal* **54** 614–626

[9] Soulaine C, Davit Y and Quintard M 2013 *Chemical Engineering Science* **96** 55–70

[10] Soulaine C, Horgue P, Franc J and Quintard M 2014 *AIChE Journal* **60** 3665–3674

[11] Davit Y, Golfier F, Latché J C and Quintard M 2019 *SIAM Journal on Applied Mathematics* **79** 1797–1822
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Comparison between DNS averaged temperatures (denoted DNS[i-1] for phase *i*) and macro-scale predictions (denoted Macro[i-1] for phase *i*) from the 3-Eq. model (Pe = 10, *t* = 0.7)</figcaption></figure>

the macro-scale PDE with various terms and complex (possibly spatially and time dependent) effective properties. Most of the time, additional terms would be specific to the pore-scale initial boundary value problem under investigation.

The proposed alternative of a N-eq description offers more flexibility. First, we can adapt the N-eq description in terms of *N*, and also the domain-decomposition topology, to characteristic length-scales of the BCs spatial distribution. Like one-equation descriptions, N-equation models can have complex mathematical structures (e.g. spatial and time dependent effective properties in convolution products). However, the choice of *N* and domain topology is already a powerful tool to approach the actual macro-scale behavior. Therefore, we tested here a first-order closure for the N-eq model which is *independent* of the actual IBVP under investigation, and gives simple expressions for the effective properties. Indeed, in the particular case used as an illustration in this paper, effective properties are matrices with constant coefficients obtained from analytical expressions (the use of symbolic programming allows to construct automatically these properties for any *N*, domain topology, and velocity field). Various tests showed that the macro-scale predictions are rather robust estimates of the actual behavior.

**References**

[1] Bear J 1972 *Dynamics of Fluids in Porous Media* (New York: Elsevier)

[2] Taylor G I 1953 *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences* **219** 186–203

[3] Koch D L and Brady J F 1987 *J. Fluid Mech.* **180** 387–403

[4] Cushman H J 1997 *The Physics of Fluids in Hierarchical Porous Media: Angstroms to Miles* (The Netherlands: Academic Publishers)

[5] Davit Y and Quintard M 2015 *Handbook of Porous Media* ed Vafai K (Taylor & Francis) pp 245–320 3rd ed

[6] Davarzani H, Marcoux M and Quintard M 2010 *International Journal of Heat and Mass Transfer* **53** 1514–1528

[7] Coats K H and Smith B D 1964 *Society of Petroleum Engineers Journal March* 73–84

[8] Mahr B and Mewes D 2008 *AIChE Journal* **54** 614–626

[9] Soulaine C, Davit Y and Quintard M 2013 *Chemical Engineering Science* **96** 55–70

[10] Soulaine C, Horgue P, Franc J and Quintard M 2014 *AIChE Journal* **60** 3665–3674

[11] Davit Y, Golfier F, Latché J C and Quintard M 2019 *SIAM Journal on Applied Mathematics* **79** 1797–1822