A (Dual) Network Model for Heat Transfer in Porous Media

Abstract

We present a dual network model to simulate coupled single-phase flow and energy transport in porous media including conditions under which local thermal equilibrium cannot be assumed. The models target applications such as the simulation of catalytic reactors, micro-fluidic experiments, or micro-cooling devices. The new technique is based on a recently developed algorithm that extracts both the pore space and the solid grain matrix of a porous medium from CT images into an interconnected network representation. We simulate coupled heat and mass transfer in these networks simultaneously, allowing naturally to model scenarios with heterogeneous temperature distributions in both void space and solid matrix. The model is compared with 3D conjugate heat transfer simulations for both conduction- and convection-dominated scenarios. It is shown to reproduce effective thermal conductivities over a wide range of fluid to solid thermal conductivity ratios with a single parameter set. Morevoer, it captures local thermal nonequilibrium effects in a micro-cooling device scenario.

1 Introduction

Heat transfer processes and non-isothermal effects are found in many technical and environmental porous-media systems. Evaporation of water from soils locally decreases the soil surface’s temperature (Shahraeeni and Or 2010), an effect also used for turbine blade cooling (Dahmen et al. 2014). Spontaneous imbibition may feature local temperature spikes on the order of several Kelvin due to surface chemistry interaction between the fluid and the solid phase (Terzis et al. 2017). Heat transfer in porous media also plays a crucial role in single-fluid-phase systems: In the face of global climate change, geothermal energy (Scheck-Wenderoth et al. 2013; Praditia et al. 2018) may become an essential contributor of clean and sustainable energy. Solar heat devices and air conditioning systems require efficient heat exchangers for which also porous media have been considered (Delalic et al. 2004; Hutter et al. 2011). The advance of micro-scale information technology throughout all areas of life comes with an increasing demand on energy for cooling. Miniaturization and integration of the cooling technology can help to increase efficiency and save cost (van Erp et al. 2020).

While energy transport in porous media (Bear and Bensabat 1989; Bear 2018) under conditions of both local thermal equilibrium and nonequilibrium is well studied on the Darcy scale (Nield and Bejan 2013), it has been considered to a far less extend in pore-network models (PNM). Often, isothermal conditions (Prat 1993; Laurindo and Prat 1996; Yiotis et al. 2001; Metzger et al. 2007) or a fixed temperature distribution (Huinink et al. 2002) are assumed. Heat transfer can be accounted for in a simplified manner by adding a “thermal network” (Figus et al. 1999) in between the network representing the void space and solving a single energy balance. Here, effective quantities are used to include the contribution of both the void and the solid phase while usually, only heat conduction is considered (Plourde and Prat 2003; Surasani et al. 2008, 2009; Médici and Allen 2013). A more sophisticated concept involves solving two separate energy balance equations for the fluid and the solid phase (Satik and Yortsos 1996). Both phases are represented by two individual networks and also advective heat fluxes within the fluid phase are considered. Coupling terms facilitate the exchange of energy between the phases. All the above-mentioned works have a common feature: they consider regular, structured porous materials with a fixed coordination number for all void pore bodies. With the rapid advance of CT-scanning (Cnudde and Boone 2013; Wildenschild and Sheppard 2013) and image processing techniques, the extraction of complex pore networks for real permeable media has considerably gained attention during the last years (Blunt et al. 2013; Xiong et al. 2016; Rabbani et al. 2020).

Recently, an efficient algorithm for the extraction of a dual network, considering both the void space and solid-matrix geometry, has been developed in Khan et al. (2019, 2020) based on a watershed algorithm (Gostick 2017). In Khan et al. (2019), the authors showed that the algorithm can reliably extract void-solid interfacial areas. They also simulated stationary diffusion separately on both void and solid networks. However, they did not perform coupled simulations considering both networks simultaneously.

In this work, we introduce a fully coupled, locally energy- and mass-conservative dual network model for the simulation of heat transfer in realistic natural porous media, allowing the consideration of pore-local thermal nonequilibrium and structural heterogeneity. Energy transfer in the dual network is modeled as a fully coupled processes using pore- and grain-local heat transfer rules derived from the analysis of local idealized but spatially resolved problems and geometrical considerations.

In contrast to embedded network techniques for simulating fractured rocks or wells (Karvounis and Jenny 2016; Koch 2020), where a network model is coupled to a homogenized model of an embedding bulk domain that may be itself a porous medium (e.g., sandstone), this work focuses on a single-scale approach. A porous medium is conceptualized as two interconnected transport systems, the void and the solid space. Both subdomains are simplified using the same model reduction technique commonly applied exclusively to the void space, where it is known as pore-network modelling (Table 1).

Table 1 List of important symbols

2 Conjugate heat transfer on the pore scale

In the following, we consider a porous domain $\Omega$ comprised of two subdomains, a fluid-filled void space domain ${\mathcal {F}}$, and a solid domain ${\mathcal {S}}$, such that $\Omega = {\mathcal {F}}\cup {\mathcal {S}}$. The systems are coupled on the void-solid interface ${\mathcal {I}}= {\mathcal {F}}\cap {\mathcal {S}}$. The composite system is assumed to be rigid and all pore spaces are filled with a fluid with viscosity $\mu$ and constant density $\rho$.

Fluid flow in the void spaces is governed by the Navier-Stokes equations. Mass and momentum balance for incompressible (single phase) fluid flow, neglecting gravity and dilatation, are given by

$$\begin{aligned} \frac{\partial (\rho {\varvec{v}})}{\partial t} + \nabla \!\cdot (\rho {\varvec{v}} {\varvec{v}}^T ) - \nabla \!\cdot ( \mu (\nabla {\varvec{v}}_{\text {f}}+ \nabla ^T{\varvec{v}}_{\text {f}}) - p_{\text {f}}{\varvec{I}} )&= 0, \end{aligned}$$

(1)

$$\begin{aligned} \nabla \!\cdot {\varvec{v}}_{\text {f}}&= 0, \end{aligned}$$

(2)

where ${\varvec{v}}_{\text {f}}$ denotes the fluid velocity, and $p_{\text {f}}$ is the fluid pressure in ${\mathcal {F}}$. The analysis and simulations in this work are restricted to the laminar flow regime.

We assume that there is no mass exchange between ${\mathcal {F}}$ and ${\mathcal {S}}$ and any fluid movement within ${\mathcal {S}}$ can be neglected (${\varvec{v}}_{\text {s}}\rightarrow 0$). This entails that if micro-pores are present in the solid, they are either assumed to be disconnected from the larger pores of the void space, or so small that fluid flow within the micro-pores can be neglected. Furthermore, we herein neglect the effect of surface roughness of the grains on flow and transport. Hence, we can prescribe ${\varvec{v}}_{\text {f}}= 0$ on the interface ${\mathcal {I}}$, corresponding to a no-flow, no-slip boundary condition for Eqs. (1) and (2).

The energy equation for an incompressible single-phase fluid in ${\mathcal {F}}$ is given by

$$\begin{aligned} \frac{\partial \left( \rho c_{\text {f}}T_{\text {f}} \right) }{\partial t} + \nabla \!\cdot \left( \rho c_{\text {f}} T_{\text {f}} {\varvec{v}}_{\text {f}} \right) - \nabla \!\cdot \left( \lambda _{\text {f}} \nabla T_{\text {f}} \right) = 0, \end{aligned}$$

(3)

with the specific heat capacity $c_{\text {f}}$, the thermal conductivity $\lambda _{\text {f}}$, and the temperature $T_{\text {f}}$ of the fluid. The energy equation in ${\mathcal {S}}$ is given by

$$\begin{aligned} \frac{\partial \left( \rho _{\text {s}}c_{\text {s}}T_{\text {s}}\right) }{\partial t} - \nabla \!\cdot \left( \lambda _{\text {s}}\nabla T_{\text {s}}\right) = 0, \end{aligned}$$

(4)

where $\lambda _{\text {s}}$ is the effective thermal conductivity and $c_{\text {s}}$ the (effective) specific heat capacity of the solid.

Finally, on the interface ${\mathcal {I}}$ between the two subdomains, suitable coupling conditions enforce the continuity of temperature and normal conductive energy fluxes,

$$\begin{aligned} T_{\text {s}}&= T_{\text {f}}, \end{aligned}$$

(5)

$$\begin{aligned} \lambda _{\text {s}}\nabla T_{\text {s}}\cdot {\varvec{n}}&= - \lambda _{\text {f}}\nabla T_{\text {f}}\cdot {\varvec{n}}, \end{aligned}$$

(6)

where ${\varvec{n}}$ is a unit normal vector on ${\mathcal {I}}$.

3 From pore network models to dual network models

Pore-network models are reduced models to simulate flow and transport processes in the pore space on the pore scale. To this end, the pore space is discretely represented as a network of pore bodies (nodes) and pore throats (links). Pore bodies are located at larger cavities in the pore space and throats are local constrictions that separate these cavities.

Starting from a three-dimensional representation of the pore space such as micro-CT image data, the network construction can be formalised by introducing a distance function that is the inverse Euclidean distance of every point in the pore space to the closest point on the void-solid interface. Pore body centers are located at the local maxima and pore throat centers are located at the saddle points. The resulting algorithm is known as the watershed algorithm (pore bodies correspond to the drainage basin and throats to the drainage divide in a watershed) (Jiang et al. 2007; Gostick 2017).

At least two more algorithms are commonly used: medial axis skeletonization (Lindquist et al. 1996; Prodanović et al. 2007) and the maximal ball algorithm (Silin and Patzek 2006; Dong and Blunt 2009). The choice of algorithm may result in rather large differences in geometry and topology of the pore network. However, resulting upscale quantities such as the intrinsic permeability of a network has been found to be relatively robust with respect to the chosen algorithm (Bhattad et al. 2011; Baychev et al. 2019). For a more in-depth discussion of these algorithms, we refer to Blunt (2017).

Recently, Khan and coworkers suggested a dual watershed algorithm that applies the network construction process analogously to the solid (Khan et al. 2019). Grains correspond to the drainage basin and the drainage divide to the contact surface between grains. Additionally, void-solid interfacial areas are extracted. The discretization process is shown in Fig. 1. The resulting dual network will also contain edges connecting pores with grains. These edges are associated with the local void-solid interfacial area and can be used to model energy and mass transfer between the two compartments.

Multi-network models have been used in many combinations to derive reduced models for flow and transport processes in the pore space (Rabbani et al. 2020). However, in the vast majority of pore-network modeling works, the solid phase is neglected entirely.

In this work, we will apply the pore network modeling idea both for the solid and the void space which are on a similar scale in many porous media. By describing the solid structure analogously to the pore structure, i.e. by using discrete approximations based on simplifying geometrical assumptions and local rules, it is possible to better incorporate the effects of the solid structure in network simulations of porous media. This choice is motivated by various interface-driven processes in porous media such as dissolution and precipitation, biomass growth, or chemical reactions, which can be more directly described by a method treating void and solid space in similar ways.

We call the resulting model a dual network model to emphasize that both the solid and the fluid phase form interconnected mass and energy transport networks of similar spatial extent but with different material properties. The transport networks exchange mass and energy at the solid-fluid interface.

4 Local thermal nonequilibrium in porous media

We consider local thermal nonequilibrium (LTNE) as a state where the fluid and solid phases in a porous medium have different average temperatures in a given control volume which contains at least one solid grain and one fluid-filled pore. LTNE in porous media occurs under various conditions. For example, in the case of processes with large heat transfer rates (chemical reactions, evaporation), interphasial heat transfer might not be efficient enough to establish a local thermal equilibrium (Oliveira and Kaviany 2001; Baytas and Pop 2002; Fichot et al. 2006; Li et al. 2007; Nuske et al. 2014; Lindner et al. 2016; Aminzadeh and Or 2017). LTNE readily occurs for transient heat transfer problems, such as hot/cold fluid injection scenarios (Rees et al. 2008), and in particular when materials with a large difference in thermal properties are involved (Truong and Zinsmeister 1978).

For the subsequent discussions and derivations, we define several dimensionless numbers that help determine whether a given process has a strong influence on the heat transfer within a porous medium. The Reynolds number describes the ratio of inertial to viscous forces

$$\begin{aligned} {\text {Re}} :=\frac{\rho v_{\text {ch}}l_{\text {ch}}}{\mu }, \end{aligned}$$

(7)

where $v_{\text {ch}}$ is a characteristic fluid velocity and $l_\text {ch}$ denotes a characteristic length scale. In this work, we consider laminar flow (moderate and low Reynolds numbers). The (heat transfer) Pèclet number describes the ratio of advective heat transport rate and conductive heat transport rate

$$\begin{aligned} {\text {Pe}} := \frac{\rho c_{\text {p},{\text {f}}} v_{\text {ch}}l_{\text {ch}}}{\lambda _{\text {f}}} = {\text {Re}}\,{\text {Pr}}, \end{aligned}$$

(8)

where ${\text {Pr}} := \mu c_{\text {p},\text {f}} \lambda _{\text {f}}^{-1}$ is the Prandtl number (${\text {Pr}} \approx 6$ for water at ${20}^{\circ }\,\hbox {C}$). Hence, even for laminar flow with moderate fluid velocities, convective heat transfer often dominates conductive heat transfer. The local Nusselt number is the ratio of convective to conductive heat transfer in the fluid at the void-solid interface. It can be formulated in terms of the ratio of the interfacial temperature gradient in normal direction and the mean cross-section-averaged channel temperature,

$$\begin{aligned} {\text {Nu}}:= -\frac{\nabla T_{\text {f}}\cdot {\varvec{n}}}{l_{\text {ch}}^{-1} ( T_{\mathcal {I}}- {\overline{T}}_{\text {f}})} = -\frac{\nabla T_{\text {f}}\cdot {\varvec{n}}}{l_{\text {ch}}^{-1} \Delta T_{\text {f}}} = \frac{{\mathfrak {h}}}{l_{\text {ch}}^{-1}\lambda _{\text {f}}}, \end{aligned}$$

(9)

where $T_{\mathcal {I}}$ is the temperature at the void-solid interface, and ${\mathfrak {h}}$ is the heat transfer coefficient. For $l_{\text {ch}}= 2R$ and fully developed laminar flow in a tube with circular cross-section (radius R) and uniform surface heat flux, ${\text {Nu}}\approx 4.36$ is a well-known analytical result (Incropera et al. 2006). However, the Nusselt number always depends on the chosen characteristic length as well as on the definition of $\Delta T_{\text {f}}$. We estimate local Nusselt numbers for a simple throat-pore-throat geometry in “Appendix 2”. Finally, the conductivity ratio is the ratio of fluid and solid effective thermal conductivities,

$$\begin{aligned} \kappa := \frac{\lambda _{\text {f}}}{\lambda _\text {s}}. \end{aligned}$$

(10)

For a water-granite system, $\kappa \approx 0.2$, while for a water-metal system, $\kappa \approx 0.002$. At the void-soil interface, void and solid normal conductive fluxes are equal, Eq. (6). Therefore, the conductivity ratio can be also used to characterize the ratio of the temperature gradients at the interface

$$\begin{aligned} \kappa = \frac{\nabla T_\text {s}\cdot {\varvec{n}}}{\nabla T_\text {f}\cdot {\varvec{n}}}. \end{aligned}$$

(11)

Hence, for $\kappa \ll 1$, spatial temperature gradients in the solid may be negligible in comparison to temperature changes in the fluid and $T_{\mathcal {I}}\approx {\overline{T}}_\text {s}$, while the opposite holds for $\kappa \gg 1$.

Lee and Vafai analyzed forced convection through homogeneous porous media on the Darcy scale (Lee and Vafai 1999). For a simple channel geometry (similar to Fig. 10) and under some assumptions regarding boundary conditions including a constant heat flux over the channel boundary, they derived an analytical solution for the temperature distribution in a porous channel in a cross-section normal to the main flow direction. They showed that LTNE in particular occurs in situations with moderate to low Biot numbers ($\text {Bi}:= {\mathfrak {h}}l_\text {ch}\lambda _\text {s}^{-1}$). Then, due to the low conductivity ratio, the temperature distribution in the solid is almost uniform. At the same time, interphasial heat transfer is the limited by the low conductivity of the fluid and cold water is readily supplied by advection (moderate to high Reynolds numbers). Therefore, the fluid temperature differs significantly from the solid temperature, particularly in the middle of channel.

5 A dual network model of conjugate heat transfer in porous media

Fig. 1

Dual network discretization of a porous medium. a Both void space and solid space are partitioned into discrete control volumes (pores in light blue and grains in light grey) which are separated from neighboring control volumes at local constrictions (pore throats in dark blue and grain contact surfaces in dark grey). An exemplary pore is marked with the letter F and an exemplary grain is marked with the letter S. b The void-solid interface I of pore F and grain S are shown in orange. Blue arrows show mass and energy transport between pores across throats, grey arrows show energy transport between grains across grain contact surfaces, and orange arrows show energy exchange between grains and pores. c An abstract dual network representation with nodes and links. Nodes (pores and grains) are associated with a volume measure and links (pore throats, solid contact surfaces, void-solid interfaces) are associated with an area measure

The discretization of void and solid space is depicted in Fig. 1. Both domains, ${\mathcal {F}}$ and ${\mathcal {S}}$ are partitioned into control volumes $K_\text {f}$ (where $\vert K_\text {f}\vert$ denotes the pore volume) and $K_\text {s}$ (where $\vert K_\text {s}\vert$ denotes the grain volume). Intersections between void and solid control volumes ${\mathbb {I}}= K_\text {f}\cap K_\text {s}$ describe the local void-solid interface with area $\vert {\mathbb {I}}\vert$. Intersections between two void control volumes ${\mathbb {T}}= K_\text {f}\cap L_\text {f}$ describe a pore throat (local constriction of the pore space) with area $A_{\mathbb {T}}$. Intersections between two solid control volumes ${\mathbb {A}}= K_\text {s}\cap L_\text {s}$ describe a contact surface (local constriction of the solid matrix) with area $A_{\mathbb {A}}$. The equations describing single-fluid-phase conjugate heat transfer introduced in Sect. 2 can be discretized based on the network structure by introducing local numerical approximations. We introduce the average pore pressure and temperature, and the average solid grain temperature as

$$\begin{aligned} p_{K_\text {f}} = \frac{1}{\vert K_\text {f}\vert }\int _{K_\text {f}}\! p_\text {f}\,\text {d}V, \quad T_{K_\text {f}} = \frac{1}{\vert K_\text {f}\vert }\int _{K_\text {f}}\! T_\text {f}\,\text {d}V, \quad T_{K_\text {s}} = \frac{1}{\vert K_\text {s}\vert }\int _{K_\text {s}}\! T_\text {s}\,\text {d}V. \end{aligned}$$

(12)

5.1 Dual network discretization of conjugate heat transfer

The fluid flow equations are discretized by integrating over a control volume $K_\text {f}$ and then simplified assuming specific local pore space geometries. Equations (1) and (2) are reduced to a local mass balance for each $K_\text {f}$,

$$\begin{aligned} \sum \limits _{{\mathbb {T}}\in {\mathcal {T}}} \rho g_{\mathbb {T}}(p_{K_\text {f}} - p_{L_\text {f}}) = 0, \end{aligned}$$

(13)

where ${\mathcal {T}}$ is the set of throats of control volume $K_\text {f}$, $g_{\mathbb {T}}$ the (advective) transmissibility coefficient for the throat ${\mathbb {T}}$, and $L_\text {f}$ is a neighboring control volume such that ${\mathbb {T}}= K_\text {f}\cap L_\text {f}$. While we neglect the influence of inertia effects in this work, non-linear flow occurring at higher Recould be included in a rather straightforward manner in future work (Balhoff and Wheeler 2009; Veyskarami et al. 2017).

Choosing from a large variety of formulations for $g_{\mathbb {T}}$ given in the literature, we consider (Patzek and Silin 2001; Valvatne and Blunt 2004)

$$\begin{aligned} g_{\mathbb {T}}= \frac{k G A_{\mathbb {T}}^2}{\mu l}, \end{aligned}$$

(14)

with $A_{\mathbb {T}}$ denoting the throat’s cross-sectional area and l is its length defined as the Euclidean distance between the two adjacent pore body centers minus their extended pore body radii, cf. (Gostick 2017). G is a dimensionless shape factor (Mason and Morrow 1991) given by

$$\begin{aligned} G = \frac{A_{\mathbb {T}}}{P_{\mathbb {T}}^2}, \end{aligned}$$

(15)

with the throat perimeter $P_{\mathbb {T}}$ , and k is a second shape factor (0.5 for circle, 0.5623 for square, and 0.6 for triangular throat cross-sections Patzek and Silin 2001; Valvatne and Blunt 2004), which we heuristically choose depending on G. The pore body flow resistance is neglected.

The fluid energy balance for the same control volume $K_\text {f}$ can be simplified to

$$\begin{aligned} \vert K_\text {f}\vert \frac{\partial \left( \rho c_\text {f}T_{K_\text {f}} \right) }{\partial t} - \sum \limits _{{\mathbb {T}}\in {\mathcal {T}}} \left\{ t_{\mathbb {T}}\left( T_{K_\text {f}} - T_{L_\text {f}}\right) - \rho c_\text {f}T^\text {up} g_{\mathbb {T}}\left( p_{K_\text {f}} - p_{L_\text {f}}\right) \right\} = Q_{K_\text {f}}, \end{aligned}$$

(16)

where $T^\text {up}$ is the temperature of cell $K_\text {f}$ or $L_\text {f}$, whichever has the largest average pressure (upwind scheme), and $Q_{K_\text {f}}$ is the sum of the energy influx at the void-solid interface, specified below. Finally, we introduced $t_{\mathbb {T}}$, the effective thermal transmissibility associated with ${\mathbb {T}}$, which is used to model microstructural effects.

For every solid control volume $K_\text {s}$, the local discrete energy balance is given by

$$\begin{aligned} \vert K_\text {s}\vert \frac{\partial (\rho _\text {s}c_\text {s}T_{K_\text {s}}) }{\partial t} - \sum \limits _{{\mathbb {A}}\in {\mathcal {A}}} t_{\mathbb {A}}(T_{K_\text {s}} - T_{L_\text {s}}) = Q_{K_\text {s}}, \end{aligned}$$

(17)

where ${\mathcal {A}}$ is the set of solid contact surfaces ${\mathbb {A}}= K_\text {s}\cap L_\text {s}$ of control volume $K_\text {s}$ with neighboring control volumes $L_\text {s}$, and $t_{\mathbb {A}}$ is the effective thermal transmissibility associated with ${\mathbb {A}}$ and models microstructural effects.

The total flux from a fluid control volume into the surrounding solid control volumes, $Q_{K_\text {f}}$ in Eq. (16), is the sum of the interface fluxes $Q_{{\mathbb {I}}}$ over all void-solid interfaces of a given pore control volume $K_\text {f}$. The total flux for a solid control volume, $Q_{K_\text {s}}$ in Eq. (17), is computed analogously for solid control volumes $K_\text {s}$. The interface flux is formulated analogously to the other local rules by

$$\begin{aligned} Q_{{\mathbb {I}}} = t_{{\mathbb {I}}}(T_{K_\text {f}} - T_{K_\text {s}}) \end{aligned}$$

(18)

with the interface transmissibility $t_{{\mathbb {I}}}$.

The (dual) network may be interpreted as representing a resistor network (shown in Fig. 2 for the example of a cubic lattice sintered sphere packing). The discretization scheme is effectively a finite volume scheme where the control volumes are pores and grains and the control volume faces comprise of throats, contact surfaces, and void-solid interfaces. Hence, the introduced transmissibility coefficients ($t_{\mathbb {T}}$, $t_{\mathbb {A}}$, $t_{\mathbb {I}}$) could be estimated by standard two-point flux approximations. This is sensible if material properties are similar in void and solid ($\kappa \approx 1$) and the fluid is at rest, since then, the microstructure has little influence on heat flow. However, for larger differences in the thermal conductivities ($\kappa \gg 1$ or $\kappa \ll 1$), the microstructure locally and nonlinearly distorts the temperature fields. The same holds true for convection-dominated heat transfer. As the microstructure cannot be resolved in the network model, microstructural information has to be incorporated into the transmissibility coefficients. In the following sections, we therefore suggest approaches for the estimation of $t_{\mathbb {T}}$, $t_{\mathbb {A}}$, and the coupling transmissibility $t_{\mathbb {I}}$ based on observations from 3D computational models.

5.2 Approximation of conductive heat fluxes

Fig. 2

Network representation of a cubic lattice packing of sintered spherical particles. Left: solid grains (grey), middle: fluid-filled pores (blue) and interfacial surface (gold), right: representation of the porous medium as a network of resistors. Blue resistors corresponds to pore-pore connections, grey resistors corresponds to grain-grain connections and orange resistors to pore-grain interfaces

To approximate effective thermal transmissibilities in the fluid and the solid, we start by investigating the heat flow in locally constricted geometries. Let us first consider the case where heat flow exclusively occurs in the void space ($\kappa \gg 1$). We introduce the cross-sectional pore area, $A_\text {f}$ (the area of the intersection of a plane through the pore center, parallel to the pore throat surface). The throat area is denoted by $A_{\mathbb {T}}= \vert {\mathbb {T}}\vert$, and $\Delta x$ is defined as the distance between pore center and throat center (or in the context of the solid, the distance between grain center and contact surface center). We can account for the local constriction of the cross-sectional area toward the throat by linear interpolation (conical shape),

$$\begin{aligned} A(x) = \left[ \sqrt{A_\text {f}} + x (\Delta x)^{-1}( \sqrt{A_{\mathbb {T}}} - \sqrt{A_\text {f}}) \right] ^2, \end{aligned}$$

(19)

where $x=0$ is the pore center and $x=\Delta x$ is the throat center. Given an area function, the effective thermal transmissibility can be estimated by using a generalized harmonic mean

$$\begin{aligned} t_{{\mathbb {T}},K_\text {f}} \approx \left[ \frac{1}{\lambda _\text {f}} \int _0^{\Delta x} \! \frac{1}{A(x)} \,\text {d}x \right] ^{-1} \mathop = \limits^{{{\text{Eq}}.(19)}} \quad \frac{\lambda _\text {f}\sqrt{A_\text {f}A_{\mathbb {T}}}}{\Delta x}. \end{aligned}$$

(20)

The effective thermal transmissibility from pore center to the neighboring pore center is found by the harmonic mean of thermal transmissibilities from each side,

$$\begin{aligned} \frac{1}{t_{\mathbb {T}}} = \frac{1}{t_{{\mathbb {T}},K_\text {f}}} + \frac{1}{t_{{\mathbb {T}},L_\text {f}}}. \end{aligned}$$

(21)

The effective thermal transmissibility for the solid contact surfaces, $t_{\mathbb {A}}$, can be computed analogously, assuming $\kappa \ll 1$. We remark that typical network extraction algorithms do not directly provide $A_\text {f}$. In this case, we estimate $A_\text {f}$ at each throat from the pore volume and the pore-throat-center distance, $A_\text {f}\approx V_\text {f}/ (2\Delta x)$. We denote with $A_\text {s}$ the corresponding area for a solid grain and with $V_\text {s}$ its volume such that the analogous estimate for the solid phase is $A_\text {s}\approx V_\text {s}/ (2\Delta x)$.

Transmissibilities estimated in this way assume that all heat flow occurs in the respective constricted geometry. However, in the presence of interphasial heat fluxes at the void-solid interface, heat flow is not exclusively occurring through the constrictions at throat and contact surfaces. In the limiting case of equal thermal conductivities in both phases ($\kappa = 1$), and when the fluid is at rest, the temperature distribution is unaffected by the geometry of internal interfaces. The heat flow across a given interface or constriction surface is proportional to the surface area projected into the flow direction. If the constriction areas are small in comparison to the interfacial area, most heat flow will occur between phases.

For a concrete example, we look at heat flow parallel to one of the coordinate axes in a cubic lattice sintered sphere packing. We can use the approach suggested in Eq. (20) to compute transmissibilities for all values of $\kappa$ by introducing effective areas ${\tilde{A}}_\text {f}$ and ${\tilde{A}}_\text {s}$,

$$\begin{aligned} t_{{\mathbb {T}},K_\text {f}} \approx \frac{\lambda _\text {f}\sqrt{{\tilde{A}}_\text {f}A_{\mathbb {T}}}}{\Delta x}, \quad t_{{\mathbb {T}},K_\text {s}} \approx \frac{\lambda _\text {s}\sqrt{{\tilde{A}}_\text {s}A_{\mathbb {A}}}}{\Delta x}. \end{aligned}$$

(22)

To estimate the effective areas, we consider heat flow through a simple geometry with a fluid and a solid phase with conductivity ratio $\kappa$ (Fig. 3).

Fig. 3

Heat conduction for different conductivity ratios $\kappa$. Temperature gradient field lines are sketched for heat flow occurring in vertical direction. Orange temperature gradient field lines cross over to the other phase. The solid matrix has a constriction on top (with contact area $A_{\mathbb {A}}$) connecting the grain to the neighboring grain. Effective areas ${\tilde{A}}$ are introduced as the area of a plane surface parallel to the contact surface that contains all temperature gradient field lines crossing the contact surface (black lines). Estimates of these effective areas for a given geometry can be used to estimate transmissibilities

We observe that ${\tilde{A}}_\text {s}\rightarrow A_\text {s}$ for $\kappa \rightarrow 0$, ${\tilde{A}}_\text {f}\approx A_{\mathbb {A}}$ for $\kappa = 1$ and that ${\tilde{A}}_\text {s}$ will converge to some constant $C_{\infty ,\text {s}}A_{\mathbb {A}}$ for $\kappa \rightarrow \infty$. The same observation holds for the solid using the inverse of $\kappa$. That means ${\tilde{A}}_\text {f}\rightarrow A_\text {f}$ for $\kappa \rightarrow \infty$, ${\tilde{A}}_\text {s}\approx A_{\mathbb {T}}$ for $\kappa = 1$, and ${\tilde{A}}_\text {f}\rightarrow C_{0,\text {f}}A_{\mathbb {T}}$ for $\kappa \rightarrow 0$. We suggest to model the effective areas over the entire range of $\kappa$ as

$$\begin{aligned} {\tilde{A}}_\text {f}&= A_{\mathbb {T}}\left( C_{\infty ,\text {f}} + \frac{ (C_{0,\text {f}} - C_{\infty ,\text {f}}) (C_{\infty ,\text {f}} - 1) }{ (C_{\infty ,\text {f}} - 1) + \kappa (1 - C_{0,\text {f}})} \right) \end{aligned}$$

(23)

$$\begin{aligned} {\tilde{A}}_\text {s}&= A_{\mathbb {A}}\left( C_{\infty ,\text {s}} + \frac{ (C_{0,\text {s}} - C_{\infty ,\text {s}}) (C_{\infty ,\text {s}} - 1) }{ (C_{\infty ,\text {s}} - 1) + \frac{1}{\kappa } (1 - C_{0,\text {s}})} \right) \end{aligned}$$

(24)

where $C_{0,\text {f}}$ and $C_{\infty ,\text {f}} \approx A_\text {f}/A_{\mathbb {T}}$ are dimensionless shape factors for the limiting cases $\kappa \rightarrow 0$ and $\kappa \rightarrow \infty$, respectively. The solid shape factors $C_{0,\text {s}}$ and $C_{\infty ,\text {s}} \approx A_\text {s}/A_{\mathbb {A}}$ correspond to the limiting cases $1/\kappa \rightarrow 0$ and $1/\kappa \rightarrow \infty$, respectively.

The conductive interface transmissibility (for resting fluid) can be estimated by

$$\begin{aligned} t_{{\mathbb {I}}} \approx \frac{ C_{\mathbb {I}}\vert {\mathbb {I}}\vert }{\Delta x_{\text {f}\text {s}}} \frac{\Delta x_{\text {f}\text {s}} \lambda _\text {f}\lambda _\text {s}}{ \Delta x_\text {s}\lambda _\text {f}+ \Delta x_\text {f}\lambda _\text {s}}, \quad \Delta x_{\text {f}\text {s}} = \Delta x_\text {f}+ \Delta x_\text {s}, \end{aligned}$$

(25)

where $x_{\text {f}\text {s}}$ is the distance from pore center to grain center and $C_{\mathbb {I}}$ is a dimensionless shape factor which is 1 for a completely flat void-solid interface and less than 1 for curved interfaces. The remaining factor is a distance-weighted harmonic mean of the thermal conductivities, where the weights account for distances travelled in the different materials.

For a known geometrical configuration, the effective thermal transmissibility estimate can be refined. As an example, we compute in “Appendix 1” the transmissibilities for the case of sintered spherical particles packed in a cubic lattice (overlapping sphere packing) by both analytical estimations and numerical simulations. Figure 4 shows the resulting thermal transmissibilities for different fractions of particle overlap computed with different methods. Moreover, we computed effective transmissibilities of cubic lattice sphere packing unit cells for different conductivity ratios (as described in “Appendix 1” and shown in Fig. 4 bottom right). It can be seen that for $\kappa \rightarrow 1$, the dependence of the transmissibility on the internal geometry decreases and vanishes in the limit ($\kappa = 1$), as expected. Furthermore, it is evident that in an intermediate regime ($0.01< \kappa < 100$) the transmissibility strongly depends on $\kappa$.

Fig. 4

Effective thermal conductivity factors for a cubic lattice packing of sintered spherical particles. Normalized solid and void thermal transmissibilities for different fraction of overlap ($\Delta x / R = 0.95$ means 5% overlap with respect to the particle radius). Top row: the solid grey line shows 1D analytical estimations (see “Appendix 1”); the dotted grey line (“cone approximation”) shows the normalized transmissibility approximated by Eq. (20) for the same geometry; the curve annotated with $V/(2\Delta x)$ estimates the pore/grain cross-sectional areas with the pore/grain volume and the pore/grain center to throat/contact center distance, $\Delta x$, instead of using the analytically computed areas; the solid black line show results obtained with 3D simulations of unit cells (see “Appendix 1”). Bottom row: the porosity of the particle packing for different overlap fractions is shown on the left; effective thermal conductivities for different conductivity ratios from 3D simulations of unit cells are shown on the right

5.3 Approximation of interfacial heat fluxes in the presence of convection

Fig. 5

Local solid heat flow patterns for the case of convective heat transfer. Heat flow (orange, lines following the largest local temperature gradient) in the solid domain for different pore throat Reynolds numbers in a rotational-symmetric wavy channel (left fluid boundary is the rotation axis). Fluid (blue) flows from bottom to top. A constant heat flux is applied on the outer solid boundary. The heat lines are the result of simulations described in detail in “Appendix 2”

Local heat flow patterns in the presence of fluid convection fundamentally differ from the case of pure conduction (resting fluid). Local temperature gradients at the void-solid interface, induced by the fluid movement, influence heat flow and completely dominate the heat flow pattern at high Reynolds numbers. To investigate heat and mass transfer in the pore space, we conducted conjugate heat transfer simulations (solving the equations in Sect. 2) in a sinusoidal wavy channel geometry for varying flow conditions. The detailed setup and results are given “Appendix 2”. The geometry and the boundary conditions for our setups are shown in Fig. 16.

When looking at the injection of a cold fluid into a heated porous medium (see Fig. 5), heat will typically flow in opposite direction of mass, from hot regions close to the heat source to cold regions close to the injection. With increasing Reynolds numbers, heat transfer at the void-solid interface increases and heat flow in the solid localizes and occurs in direction of the local interface. This effect is strongest where local Reynolds numbers are highest, that is in the constriction of the fluid pore space (throats, see Fig. 10). By simulation of mass and energy transport in wavy fluid channels, Wang and Vanka (1995) show that local Nusselt numbers in constrictions can easily exceed those in the cavities by an order of magnitude. This motivates a throat-centered approach to modeling convective heat transfer in the network model. Moreover, when looking at heat flow in the solid (Fig. 5), there are certain configurations where the average temperature difference between a solid grain and a fluid pore in downstream direction does not characterize the interfacial heat flux well. In fact, for lower Reynolds numbers (right image in Fig. 5) these temperatures may be equal although the interface flux is not, a configuration which we observed in our numerical experiments in wavy channels (“Appendix 2”). On the other hand, the temperature difference to an upstream pore may overestimate the average solid-fluid temperature difference along the void-solid interface.

We therefore propose to estimate effective heat transfer coefficients based on the average throat temperature defined as the arithmetic mean of the adjacent pore temperatures $T_{\mathbb {T}}= 0.5(T_{K_\text {f}} + T_{L_\text {f}})$. We define a throat is in thermal contact with a solid grain if both adjacent pores have a non-zero interfacial area with the solid grain. With each throat, we associate an upstream and a downstream void-solid interfacial area $A_\text {up}$ and $A_\text {dn}$. The upstream area is computed as the total interfacial area of the upstream pore $K_\text {f}$ with a given solid grain, divided by the number of throats which are both connected to the same upstream pore and in thermal contact with the same solid grain. The downstream interfacial area is defined analogously. Furthermore, we denote with $d_{\mathbb {T}}$, the distance between the throat center and the center of the coupled solid grain.

Let us assume that the interfacial surface between throat and solid grain is given by $I_{\mathbb {T}}$ and its area is $\vert I_{\mathbb {T}}\vert = A_\text {up} + A_\text {dn}$ and $T_{\mathbb {I}}$ the average interface temperature. Then, we can formulate a discrete heat flux between throat and grain by introducing a local Nusselt number $\text {Nu}_{\mathbb {T}}$ such that

$$\begin{aligned} Q_{{\mathbb {T}}\rightarrow I} = \int _{I_{\mathbb {T}}} \lambda _\text {f}\nabla T_\text {f}\cdot {\varvec{n}} \text {d}A \approx \frac{\vert I_{\mathbb {T}}\vert \lambda _\text {f}}{d_{\mathbb {T}}} (T_{\mathbb {I}}- T_{\mathbb {T}})\text {Nu}_{\mathbb {T}}, \end{aligned}$$

(26)

where ${\varvec{n}}$ is a unit normal vector on $I_{\mathbb {T}}$ directed toward ${\mathcal {S}}$, and $T_{\mathbb {I}}$ is the average temperature on the interface. Likewise, we can formulate a discrete heat flux on the solid side of the interface by introducing a local solid thermal conductivity factor $b_{\text {s},{\mathbb {T}}}$ such that

$$\begin{aligned} Q_{\text {s}\rightarrow I} = -\int _{I_{\mathbb {T}}} \lambda _\text {s}\nabla T_\text {s}\cdot {\varvec{n}} \text {d}A \approx -\frac{\vert I_{\mathbb {T}}\vert \lambda _\text {s}}{d_{\mathbb {T}}} (T_{\mathbb {I}}- T_{K_\text {s}})b_{\text {s},{\mathbb {T}}}, \end{aligned}$$

(27)

The factor $b_{\text {s},{\mathbb {T}}}$ describes how changes in the local three-dimensional heat flow map relate to the temperature difference in the solid $(T_{\mathbb {I}}- T_{K_\text {s}})$. Due to the energy flux continuity, Eq. (6), these fluxes have to be equal and the average interface temperature $T_{\mathbb {I}}$ can be eliminated:

$$\begin{aligned} Q_{{\mathbb {T}}\rightarrow \text {s}} \approx t_{\mathbb {I}}(T_{\mathbb {T}}- T_{K_\text {s}}), \quad t_{\mathbb {I}}= \frac{\vert I_{\mathbb {T}}\vert }{\Delta d_{\mathbb {T}}} \frac{\lambda _\text {f}b_{\text {s},{\mathbb {T}}}}{\kappa + b_{\text {s},{\mathbb {T}}} \text {Nu}_{\mathbb {T}}^{-1}} = \frac{\vert I_{\mathbb {T}}\vert }{\Delta d_{\mathbb {T}}} \lambda _{\text {conv}}, \end{aligned}$$

(28)

where we introduce $\lambda _{\text {conv}}$ as the effective “conductivity” for the local convective heat transfer (units of $\hbox {W m}^{-1} \hbox { K}^{-1}$). The interfacial heat flux in terms of the average solid and fluid temperature involves estimates for the interfacial temperature gradients from both sides of the interface. The flux computation may be simplified depending on the dimensionless quantities $\kappa$, $\text {Nu}_{\mathbb {T}}$, $b_{\text {s},{\mathbb {T}}}$. For instance, for the common case that the solid thermal conductivity is much larger than the fluid thermal conductivity ($\kappa \ll 1$), the expression may be simplified to $\lambda_{{\text {conv}}} \approx {\lambda_{\text {f}}}{\text {Nu}}_{\mathbb {T}}$.

When implementing Eq. (28) for a dual network model, we experienced that the scheme results in spurious temperature over- and undershoots. These instabilities are yet to be closely investigated. More importantly, distributing the interfacial heat transfer equally between upstream and downstream interface would commonly overestimate the flux from solid grain to downstream fluid pore and underestimate the flux to the upstream pore due to the local heat flow map shown in Fig. 5. We therefore propose a split of the heat flux by using the definition of the throat temperature

$$\begin{aligned} Q_{{\mathbb {T}}\rightarrow \text {s}} \approx t_{\mathbb {I}}(T_{\mathbb {T}}- T_{K_\text {s}}) = \frac{t_{\mathbb {I}}}{2}(T_{K_\text {f}} - T_{K_\text {s}}) + \frac{t_{\mathbb {I}}}{2}(T_{L_\text {f}} - T_{K_\text {s}}), \end{aligned}$$

(29)

where the first term is added as source term to the energy balance of control volume $K_\text {s}$, and the second term to the energy balance of $L_\text {s}$. Since the downstream pore is expected to be slightly warmer, this split redistributes the source in favor of the upstream control volume.

Local Nusselt numbers and solid conductivity factors generally depend on the flow conditions, the geometry of the pore space, and thermal properties of fluid and solid. We estimated Nusselt numbers and effective conductivities in a sinusoidal wavy channel for varying flow conditions, heat capacity ratios, and geometries. The detailed setup and results are given “Appendix 2”. Figure 6 shows streamline visualizations for several Reynolds numbers. For $\text {Re}_{\mathbb {T}}> 10$ we observed laminar recirculation zones caused by inertial effects. These recirculation zones enhance convective heat transfer. Due to the geometry of the pore space, larger Reynolds numbers may locally occur in throats even when the average sample velocity is small. The convective transmissibility is strongly dependent on the local Reynolds number. Motivated by the numerical results in the wavy channel (Fig. 17), and for the purpose of the proof of concept in this paper, we propose a simplified law for $\lambda _{\text {conv}}$ depending on the local Reynolds number in the throat $\text {Re}_{\mathbb {T}}$ and a fitting factor $\epsilon$,

$$\begin{aligned} \lambda _{\text {conv}} \approx \epsilon \text {Re}^{0.4}, \end{aligned}$$

(30)

where $\epsilon \approx 0.2$, 0.75, and 0.9 for the conductivity ratios $\kappa = 0.26$, 0.033 and 0.0033, respectively. We also observed a dependence on the pore geometry (in Fig. 18 we varied the pore radius at constant throat radius), in particular in the region of moderate Reynolds numbers where the heat flow profile significantly changes orientation, due to an increasingly more efficient heat transfer. For simplicity, we chose to use a fit that matches the average throat and pore radii of the Berea sandstone sample used in the numerical results section.

Fig. 6

Streamlines for different $\text {Re}_{\mathbb {T}}$ in an axisymmetric channel. The Reynolds number $\text {Re}_{\mathbb {T}}$ is defined with mean velocity and diameter in the throat. Boundary conditions are shown in Fig. 16. Color shows temperature and is rescaled to the minimum and maximum temperature for each case to increase contrast. Corresponding Péclet numbers from bottom to top are $\text {Pe}_{\mathbb {T}}= (12.4, 124, 620)$ and $\kappa = 0.003$. In all cases $\text {Pr}= 6.2$. For the low Reynolds number case, heat conduction has a strong influence and leads to local thermal equilibrium after only one pore while for higher $\text {Re}_{\mathbb {T}}$, we observe local thermal nonequilibrium over longer distances. For higher $\text {Re}_{\mathbb {T}}$, laminar recirculation zones appear in the pore cavities—a process driven by inertial forces

5.4 Implementation aspects of the dual network model

The dual network model has been implemented in the open-source numerical software framework DuMu^x (Koch et al. 2020) which already contains a variety of pore-network model implementations in a development version (Weishaupt et al. 2019; Weishaupt 2020). Moreover, it provides tools to implement multi-domain simulations in which coupled equations in multiple domains can be solved with possibly different computational meshes or discretization schemes. Dual networks, extracted from binary images using the open-source image analysis tool PoreSpy (Gostick et al. 2019; Khan et al. 2019), are split in DuMu^x into two separate meshes for void and solid domain (using the Dune grid interface (Bastian et al. 2008) implementations dune-subgrid (Gräser and Sander 2009) and dune-foamgrid Sander et al. 2017). The coupling information at void-solid interfaces is extracted and allows to reconnect the meshes and implement the coupling conditions. While this weak coupling approach in the software implementation is not strictly necessary for the presented model and results in more verbose code and computational overhead with respect to a minimal implementation, it is more flexible and allows to exchange the networks models individually. For example, this can be used to add additional physical processes (such as multi-phase, multi-component flow in the void domain) in the future. The discrete equations are assembled in residual form

$$\begin{aligned} \begin{bmatrix} J_\text {f}&{} C_{\text {f}\rightarrow \text {s}} \\ C_{\text {s}\rightarrow \text {f}} &{} J_\text {s}\end{bmatrix} \begin{bmatrix} \Delta u_\text {f}\\ \Delta u_\text {s}\end{bmatrix} = \begin{bmatrix} r_\text {f}\\ r_\text {s}\end{bmatrix}, \end{aligned}$$

(31)

where $J_\text {f}$ and $J_\text {s}$ are the Jacobians of the discrete PDE systems (see Sect. 5) for subdomain ${\mathcal {F}}$ and ${\mathcal {S}}$, respectively. The coupling conditions are assembled in form of the coupling Jacobian block $C_{\text {f}\rightarrow \text {s}}$ containing the derivatives of residuals of one domain with respect to degrees of freedom of the other domain, $C_{\text {f}\rightarrow \text {s}} = \frac{\partial r_\text {f}}{\partial u_\text {s}}$. The unknown vector $u_\text {f}$ contains pore fluid pressures and temperatures and $u_\text {s}$ contains the solid grain temperatures. The symbols $\Delta u_\text {f}$ and $\Delta u_\text {s}$ denote unknown-vector updates with respect to some intial guess $[ u_{\text {f},0}, u_{\text {s},0} ]^T$. The linear system in Eq. (31) is solved with a direct solver. The convection system with its nonlinear dependence of the convective transmissibility on the local Reynolds number results in a nonlinear system. Here, we linearize, solve Eq. (31) and update the initial guess in every iteration of a Newton method.

6 Numerical results and discussion

In this section, the novel dual network is applied for two types of problems: heat conduction and forced convective heat transfer in porous media. First, we consider heat conduction in a cubic lattice packing of sintered spheres and in a sample of Berea sandstone. Network computations are compared with three-dimensional simulations of heat conduction. Secondly, we consider the problem of forced convection of a cooling fluid through a porous medium heated from below. Again, we compare network computations with 3D resolved simulations.

6.1 Heat conduction in porous media

6.1.1 Sintered sphere packing

A periodic cubic lattice packing of sintered spheres allows to perform a detailed analysis of the novel model’s accuracy as all geometrical parameters and quantities can be determined analytically. Three-dimensional simulations of heat conduction in a unit cell (containing 1/8 of a sintered sphere) allow us to separately compute conductive fluxes through fluid throats, solid contact areas, and over the void-solid interface. The derived heat fluxes together with analytical geometrical parameters (distances and areas) and the effective conductivity model proposed in Sect. 5.2 are used to parameterized a dual network model of the geometry shown in Fig. 2. The results of this parameter fitting are shown in “Appendix 1”. For comparison, we also used the open-source software PoreSpy (Gostick 2017; Gostick et al. 2019; Khan et al. 2019) to extract a dual network from artificially created 3D voxel images of the sintered sphere packing (image resolution such that the sphere radius $R = 50~\text {px}$). Comparison of the dual network parameterized with analytical distances and areas and the network parameterised with algorithmically extracted quantities helps to verify correctness of the extracted microstructural information. Simulation results for effective thermal conductivities of sintered sphere packings with different amounts of overlap are shown in Fig. 7. The dual network with exact areas and distances matches the 3D simulation results well over a wide range of $\kappa$. Deviations result from errors in the conductivity model (Sect. 5.2). The results for the network parameterised with by areas and distances extracted from binary images are very similar, underlying the high accuracy of the extraction algorithm (Khan et al. 2019). The small deviations can be attributed to the discrete nature of the voxel images serving as input the extraction process. We note that these results could only be obtained by enabling PoreSpy ’s marching-cube algorithm which considerably increases the accuracy of the extracted pore-throat and interconnection areas at the cost of increased run time for the extraction process.

Fig. 7

Effective thermal conductivity of sintered sphere packing. Effective thermal conductivities of a lattice cube sintered sphere packing. Results from 3D unit cell simulations (“Appendix 1”) and dual network simulations. Transmissibility parameters for void, solid and void-solid connections for the dual network are obtained by fitting the proposed conductivity model (Sect. 5.2) to 3D simulation data on unit cells (“Appendix 1”). The employed parameters are $C_{\infty ,\text {f}} = A_{\mathbb {T}}A_\text {f}^{-1}$, $C_{0,\text {f}} = 0.1$, $C_{\infty ,\text {s}} = 1.45 A_{\mathbb {A}}A_\text {s}^{-1}$, $C_{0,\text {f}} = 0.45$, $C_{\mathbb {I}}= 0.9$. While for the case “analytic”, we computed areas and distances for the transmissibilities analytically, the case “extraction algorithm” uses areas and distances as extracted by PoreSpy from binary images of sphere packings

6.1.2 Berea sandstone sample

For the second test case, we have used the CT-image data of a Berea sandstone sample¹ with a resolution of $5.345 \times 10^{-6}\hbox {m}$ provided by Dong and Blunt (2009) . We first cut out a cubic domain of $200^3$ voxels with a side length of $L^{\mathrm {Berea}} = 1.069\times 10^{-1} \hbox {m}$ and a porosity of 0.197 from the center of the image and used PoreSpy for extracting a dual network with 556 void and 560 solid pores, respectively. The 3D voxel image and the resulting dual network are shown in Fig. 8.

Fig. 8

Extracted dual network from the Berea sandstone sample. The left figure shows the CT-scan voxel image of the solid phase in grey and the corresponding fluid-phase pore network embedded into its void space in blue. The voxel image is cubic but a small part has been cut away here for better visualisation. The center image shows the three-dimensional void space of the Berea sample and the corresponding solid-phase pore network. The image on the right depicts the dual pore network with fluid pores in blue and solid pores in grey. Interconnecting throats are shown in purple

Fig. 9

Effective thermal conductivity of Berea sandstone sample. Left: effective thermal conductivity of the sample for different conductivity ratios in three coordinate directions. The thermal conductivity is found to be isotropic. Comparison of data obtained by 3D simulations on the original binary image, and result of the dual network approach on a network extracted by PoreSpy. Individual conductivities are parameterised by the proposed model with $C_{\infty ,\text {f}} = A_{\mathbb {T}}A_\text {f}^{-1}$, $C_{0,\text {f}} = 0.1$, $C_{\infty ,\text {s}} = 0.5 A_{\mathbb {A}}A_\text {s}^{-1}$, $C_{0,\text {f}} = 0.4$, $C_{\mathbb {I}}= 0.52$. Right: heat flow streamlines from 3D simulation through a thin slab of the sandstone sample for $\kappa = 1\times 10^{-6}$. White circles mark contact regions where heat flow is evidently forced through small constrictions reducing the effective (apparent) thermal conductivity of the porous medium

In order to asses the accuracy of Eqs. (22) and (25) for heat conduction in a natural porous medium in absence of fluid flow (${\varvec{v}}_\text {f}= 0$), we first created numerical reference solutions by solving Eqs. (3) and (4) on a structured 3D Cartesian grid with uniform, axis-aligned hexagonal grid cells coinciding with the CT-image’s voxels. A thermal conductivity of $\lambda _f$ or $\lambda _s$ was assigned to grid cells corresponding to void or solid voxels, respectively. A two-point-flux cell-centered finite volume scheme implemented in DuMu^x was used to create the steady state-solution. We subsequently applied a temperature gradient in x, y and z-direction by setting a fixed temperature at the corresponding inlet and the outlet boundaries and using $\nabla T_\text {s}\cdot {\mathbf {n}} = 0$ or ${\varvec{v}}_\text {f}= 0$ at the remaining lateral boundaries. After each simulation run, the effective macro-scale thermal conductivity $\lambda ^{\mathrm {eff}}_i$, $i \in \{x,y,z\}$, was determined by integrating the total heat flow leaving the domain, $Q^{\mathrm {tot}}_i$, and evaluating $\lambda ^{\mathrm {eff}}_i = - Q^\mathrm {tot}_i (A^{\mathrm {Berea}} \nabla T_i )^{-1}$, where $A^{\mathrm {Berea}} = 1.143 \times 10^{-6}\hbox {m}^{2}$ is the area over which the heat flux leaves the domain and $\nabla T_i$ is the applied temperature gradient. This procedure was repeated for a range of different conductivity ratios. We then used the dual network model to estimate $\lambda ^{\mathrm {eff}}_i$ based on the same principle as described above. The individual transmissibilities are estimated as described in Sect. 5.2 and the parameter were tuned manually to obtain a good fit over the entire range of thermal conductivity ratios. The same parameters were used for all three coordinate directions and areas and distances to compute the transmissibilities are those extracted by PoreSpy. The results are shown in Fig. 9. An excellent agreement over the entire range of $\kappa$ could be achieved and the sample shows an isotropic effective conductivity in both the reference simulation and with the dual network model. The subtle trend that the effective conductivity in z-direction seems to be slightly higher than the effective conductivity x-direction for low $\kappa$ is also recovered with the dual network approach.

6.2 Convective heat transfer with external heat flux

We simulate forced convection through the Berea sandstone sample acting as a micro-cooling device and compare the novel network model’s results with reference solutions obtained by three-dimensionally resolved pore-scale conjugate heat transfer simulations performed with the open-source CFD toolbox OpenFOAM (Jasak 2009). All simulations are either run until a steady-state is established (OpenFOAM), or we solve for the steady-state directly (network model).

Fig. 10

Numerical scenario: a porous micro-cooling device. Left, domain and boundary conditions. Top, front and back are insulated (no-flow, no-slip boundary conditions). The bottom plate has a constant temperature of ${400}\hbox { K}$, and cold fluid at ${300}\hbox { K}$ enters on the left boundary and leaves the domain on the right. The flow is driven by the pressure gradient imposed by the left and right boundary conditions. Right, 3D simulation result with OpenFOAM: streamline visualisation of fluid flow through the Berea sandstone sample for a prescribed pressure gradient of ${9.35 \times 10^{6}} \hbox {Pa m}^{-1}$. Color shows velocity magnitude normalised by the mean velocity magnitude. Local velocity peaks in throats can exceed the mean velocity by one or more orders of magnitude, cf. (Muljadi et al. 2016)

The scenario is sketched in Fig. 10. The Berea sandstone cube acts as a heat exchanger. It is placed onto a hot bottom surface with constant temperature and cold fluid flow is forced from left to right. We note that the scenario is deliberately kept simple to ease model comparison and the geometries and parameters chosen make no attempt in representing a viable candidate for a particularly efficient micro-cooling device.

Since we focus on convective heat transfer, the hydraulic transport properties of the porous medium have a crucial influence on the overall system behavior. In a preliminary study, we compared permeabilities computed with the network model to published data obtained with 3D Lattice-Boltzmann simulations (Dong and Blunt 2009) (see supplementary material S1). In z-direction of the Berea sandstone sample, we obtained a very good agreement ($K = {1.18 \times 10^{-12}}\,\hbox {m}^{2}$). Therefore, we choose this flow direction for the subsequent numerical analysis of forced convection.

The fluid’s parameters are given by $\mu = {1\times 10^{-3}}\hbox {Pa s}$, $\rho = {1000}\,\hbox {kg m}^{-3}$, $\lambda _\text {f}= {0.679}\,{\hbox {W K}^{-1}\hbox {m}^{-1}}$, $c_\text {f}= {4.2}\,\hbox {kJ kg}^{-1}\hbox { K}^{-1}$, resulting in a Prandtl number of $\text {Pr}= 6.19$. Constant fluid parameters—although temperature ranges are high—are chosen not to add further complexity to the model comparison. We run simulations for three different solids: solid A (aluminium-like, high conductivity) with $\lambda _\text {s}= {205}\,\hbox {W K}^{-1}\hbox {m}^{-1}$, solid B (lead-like, intermediate conductivity) with $\lambda _\text {s}= {26}\,\hbox {W K}^{-1}\hbox {m}^{-1}$ and solid C (granite-like, low conductivity) with $\lambda _\text {s}= {2.6}\,\hbox {W K}^{-1}\hbox {m}^{-1}$, resulting in $\kappa \approx 0.2612$, $\kappa \approx 0.033$, $\kappa \approx 0.0033$, respectively. The solid heat capacity and density do not affect the simulation results since we simulate steady-state conditions.

We performed 3D pore-scale conjugate heat transfer simulations of the described scenario using the ChtMultiRegionFoam module of OpenFOAM to obtain a numerical reference solution. Due to the small extent of the domain and the difference in numerical schemes care has to be taken to implement matching boundary conditions. Details on the chosen grid and the implementation of the boundary conditions is described in detail in “Appendix 3”.

6.2.1 Numerical results for the micro-cooling device scenario

Fig. 11

Simplified sketch of the micro-cooling device scenario with global energy fluxes

For evaluating the results and for the model intercomparison, we compute the energy fluxes $Q^h_{\text {in,ad}}$, $Q^h_{\text {out,ad}}$, $Q^h_{\text {out,cond}}$, $Q^h_{\text {in,cond}}$, and $Q^h_{{\mathbb {I}}}$ as sketched in Fig. 11 and explained in more detail in “Appendix 4”. Energy enters the domain by advection on the inlet and conduction at the bottom plate. It leaves the domain by advection on the outlet and conduction at the inlet (since the external fluid temperature $T_{\mathrm {in}} = {300}\hbox { K}$ is lower than the internal fluid and solid temperatures raised by the heating from below). The inlet flux at the bottom plate is forced to penetrate over the void-solid interface. Finally, we computed the total mass flow rate at the inlet, ${\dot{m}}$. The results for both 3D simulation and dual network simulations are shown in Table 2. The different cases differ in the thermal conductivity ratio $\kappa$ (we varied the solid thermal conductivity), and in the applied pressure gradient across the sample geometry.

Prior to discussing more details, we need to address some possible limitations of our approach followed here: we performed all our 3D simulations on a grid featuring $200^3$ cells (corresponding to the original voxel image) without any additional refinement for reasons of limited computational resources. It this therefore unlikely that full grid convergence is already reached. However, we briefly assessed the influence of a single step of uniform grid refinement on a $100^3$ voxel center-cutout of the domain and found a decrease of less than ${7}{\%}$ for all balanced energy fluxes introduced above and the total mass flow rate. This goes well in line with the results of Muljadi and coworkers who report a change in predicted permeability of around $-{5}{\%}$ for a bead pack sample after the grid had been refined once (starting from the CT-image’s original resolution) (Muljadi et al. 2016). We therefore expect a slight overestimation of the pore-local and global mass flow rates for our 3D simulation results. Since local LTNE becomes especially relevant for higher Reynolds numbers, we tested a range of different mass flow rates, keeping in mind that we might already go beyond the realm of Darcy-type creeping flow (Muljadi et al. 2016). In contrast to the 3D model, our network model currently does not consider inertia effects and therefore the mass flow rate scales linearly with the applied pressure gradient.

We present the results for two different pressure gradients and the three different materials, A-C, as defined above in Table 2: Cases A1-C1 feature $\partial _z p = {9.35 \times 10^{6}}\hbox {Pa}$ while $\partial _z p = {9.35\times 10^{5}}\hbox {Pa}$ holds for cases A2-C2. Last, case A3 deals with an increase pressure gradient $\partial _z p = {4.68\times 10^7}\hbox {Pa}$ and the metal-like solid material.

We observe two general trends: the heat exchange rate $Q^h_{{\mathbb {I}}}$ increases non-linearly with the mass flow rate as seen from the OpenFOAM simulations. This behaviour is also captured by the network model. For a given mass flow rate, an increase in solid thermal conductivities leads to increased heat exchange rates, too.

For cases A1-C1, we see that the percentage deviation $\Delta \%$ between the network model’s results and the ones obtained with OpenFOAM for all measures defined above is less than ${10}{\%}$ except for case C1 where $Q^h_{\text {out,cond}}$ deviates noticeably. This could be related to the chosen type of Robin boundary condition which generally seems to underestimate the conductive boundary energy flux. The pore-local boundary areas were only estimated based on $R_p$ and we did not perform any further fitting here. In general, we can observe higher deviations for lower exchange rates. This behaviour is also related to the approximate (and thus, potentially inaccurate) nature of the network model’s boundary conditions, the numerical values of which become close to the internal heat exchange rate itself.

Cases A2-A3 in Table 2 feature a pressure gradient in z-direction of $\partial _z p = {9.35\times 10^5}\hbox {Pa}$ leading to maximal throat-local Reynolds numbers $\text {Re}_{\mathbb {T}}$ in the network model of up to 3.58 while the corresponding mean value is 0.24. This means that the assumption of creeping, Poiseuille-type flow in the pore throats should be generally met. Comparing the mass flow rates ${\dot{m}}$ of OpenFOAM and the network model, we see that the latter underestimates the sandstone sample’s permeability by around ${18}{\%}$, resulting in a comparable deviation for the convective heat flow rates $Q^h_{\text {in,ad}}$ and $Q^h_{\text {out,ad}}$. Returning to mass flow rates, we see that the 3D simulation exhibits a sublinear scaling of ${\dot{m}}$ with respect to $\partial _z p$ for cases A1-A3 and A3 due to the onset of non-Darcy flow at increased Reynolds numbers (Muljadi et al. 2016). his behaviour cannot be reflected by the network model due to the linear nature of the Poiseuille-type law for flow resistance. We therefore observe a seemingly close fit of ${\dot{m}}$ for cases A1-A2 where the network model’s general underestimation of the flow is outbalanced by its disability to reflect sublinear mass flow rate scaling. This trend becomes even more clear for case A3 where the network clearly over-predicts the sample’s permeability.

To clarify if the underestimated mass flow is relevant for the deviation in exchange rates, we linearly scaled the global pressure gradient of the network model such that ${\dot{m}}$ is comparable to the one of OpenFOAM in cases B2 and A3, marked with a $(*)$ in Table 2. However, the adjustment of the mass flow rate only slightly improves the results of the interfacial flux $Q^h_{{\mathbb {I}}}$ for case B2 because the conductive boundary fluxes have a relatively high influence and are not affected by the increase of ${\dot{m}}$. The adjusted mass flow rate even leads to higher deviations for case A3. This shows that a simple global scaling of the mass flow rate is not sufficient. Inertial effects act at a highly local scope, influencing both the flow resistance and the heat exchange for higher Reynolds numbers. Of course, this needs further attention and network models including inertia terms such as described in Balhoff and Wheeler (2009), Veyskarami et al. (2017) should be considered in future work.

Fig. 12

Vertical temperature distribution in micro-cooler scenario. Temperature distribution for case A3 ($\kappa = {3.3\times 10^{-3}}, \partial _z p = {4.675\times 10^7}{\hbox {Pa m}^{-1}}$) for both 3D simulations and the dual network model (NM). Mean temperatures in z-y-slabs over the x-direction (height above the hot bottom plate). The scatter plot shows the discrete temperature values in pores (blue) and grains (orange) from the dual network simulation

Fig. 13

3D temperature map in micro-cooler scenario. Temperature map for case A3 ($\kappa = {3.3\times 10^{-3}}, \partial _z p = {4.675\times 10^7}\hbox {Pa m}^{-1}$). The visualisation shows an overlay of the 3D simulation result with OpenFOAM with the dual network simulation result. The fluid domain is shown on the left, the solid domain on the right. Pores and grains of the dual network model are visualised as spheres and their size is proportional to the respective pore or grain volume

Table 2 Energy fluxes for cooling scenario

In Fig. 12, the vertical temperature distribution in both the fluid and the solid phase is given. Both temperatures decrease with increasing vertical distance from the heater plate. LTNE of up to ${10}\hbox { K}$ can observed and the network model seems to slightly overestimate the mean temperature difference between solid and fluid.

Almost all solid pores fall into a narrow temperature band while there is more variation within the fluid pore temperatures. This could be explained by the high solid conductivity smoothing out temperature gradients such that the mean solid temperature varies smoothly over height. In contrast, the fluid mean temperature shoes fluctuations, probably caused by the rather strong heterogeneity of the local fluid flow field, affected by the sample’s tortuosity and the presence of preferential flow paths. In Fig. 13, we visualise the 3D solid and fluid simulation results, and their respective network representations, separately. There is a close visual match, with some local deviations in the fluid phase temperatures. This demonstrates the network model’s ability to resolve the three-dimensional temperature field.

Finally, we note that the average computation time for one 3D calculation in Table 2 was 187.5 core hours (parallelized on 125 cores²). The average computation time for one dual network computation was several seconds on a regular laptop (single core execution) without optimizing the implementation for speed. This corresponds to a more than 100 000-fold reduction in CPU time.

7 Summary and conclusion

We presented a dual network approach to simulate heat transfer in porous media under conditions of local thermal nonequilibrium at the microstructural level. The approach allows to directly consider effects of structural features of the solid matrix, such as interphasial contact areas, local heterogeneities, and connectivity. Our monolithic coupling strategy inherently guarantees the local conservation of mass and energy. The model concept opens the way for simulations incorporating additional physical processes at the micro-structural level. This is of great interest for the understanding of systems including chemical reactions, combustion processes, evaporation, or dissolution and precipitation processes in porous media, in particular where these effects occur in small-scale systems or at material interfaces or surfaces. As for other applications of pore-network modelling, the quality of model predictions hinges on the accuracy of the derived local rules for heat and mass transfer, and is dependent on how well the studied porous medium matches the assumptions put in the derivations of these rules. We suggested local rules based on three-dimensionally resolved simulations of pore-scale flow and heat transfer.

We performed two types of numerical experiments: firstly, pure conduction with a fluid was examined for perfectly regular porous media (cubic lattice sintered sphere packing) and a sample of Berea sandstone. We found a very good agreement between the dual network’s results and 3D resolved numerical reference solutions for both types of permeable media for a wide range of fluid and solid thermal conductivity ratios. We want to point out that the analysis of porous media samples in terms of the proposed conductivity model parameters (Sect. 5.2) may pose a new way to characterize and distinguish thermal properties of samples in addition to other macro-scale indicators such as permeability and porosity.

Secondly, we considered forced heat convection with the Berea sandstone material in the role of a heat exchanger device. We performed detailed 3D conjugate-heat-transfer simulations using OpenFOAM for comparison with the dual network’s results, considering a range of different mass flow rates and material properties. The novel network model is able to reproduce the 3D results and captures changes with respect to different Reynolds number and solid materials while being computationally much more efficient. While it only takes seconds to run the network model on a regular laptop computer, the wall clock time for a forced convection OpenFOAM simulation was more than an hour on 125 cores on a cluster. In future work, we aim at integrating inertia effects into our model concept to make it applicable over a larger range of Reynolds numbers.

We believe that the presented model for laminar single-fluid-phase flow and heat transfer in porous media is a good starting point for the development of models including multi-phase multi-component transport and an important step toward more efficient model concepts for coupled systems including microstructural processes. To this end, the effect of a two-phase system within a pore on heat conduction, and the effect of structural changes due to precipitation and dissolution, thermal expansion, or grain swelling need to be investigated.

Such models may have the potential to provide additional insight when modeling fuel cell diffusion layers, porous catalysts, micro-cooling devices, evaporative cooling devices, and soil-atmosphere interaction.

Appendices

Appendix 1

1.1 Effective conductive thermal transmissibility for sintered spherical particles

We estimate the effective thermal transmissibilities $t_{\mathbb {T}}$ and $t_{\mathbb {A}}$ for a (sintered) regular cubic lattice packing of spherical particles. The spheres overlap, thus creating circular contact surfaces, ${\mathbb {A}}$. The sphere radius is denoted by R and $\Delta x$ denotes the distance from sphere center to the center of ${\mathbb {A}}$. We computed both analytical and numerical approximations.

Analytical transmissibility models We estimate $t_{\mathbb {A}}$ and $t_{\mathbb {T}}$ by one-dimensional approximations of the steady-state heat equation on the line from pore to throat centers, or grain to contact surface centers, respectively. The transmissibilities $t_{\mathbb {A}}$ and $t_{\mathbb {T}}$ can be estimated by integration of the inverse area (generalised harmonic mean). We note that this approach neglects radial heat transport and therefore only provides meaningful approximations if the rate of change of the cross-sectional area in flow direction is small. The procedure can be simplified by interpolating the area’s circle-equivalent radius linearly in the between pore center and throat center. This approach has been described in Sect. 5.2 and does not require knowledge of the exact geometry.

Fig. 14

Cubic lattice packing of sintered spherical particles. Grain (left) and pore (right) cross-sectional area shapes along the principal lattice direction

The cross-sectional areas for the cubic lattice sintered sphere packing have been computed as follows. For the solid, the cross-section starts as a disk with cut circle segments at the particle center ($x=0$, Fig. 14S), becomes eventually circular and shrinks until it reaches the circular contact surface area at $x=\Delta x$ (A in Fig. 14S):

(32)

where $x_c = \sqrt{R^2 - (\Delta x)^2}$, $r = \sqrt{(R - x)(R + x)}$, and

(33)

is the area of the overlapping circular segment (red shape in Figs. 14S and 14F). The corresponding effective thermal transmissibility from particle center to contact center is estimated by integration (generalised harmonic mean),

$$\begin{aligned} t_{{\mathbb {A}},K_\text {s}} = \lambda _\text {s}\left[ \int _0^{x_c}\! \frac{1}{A_\text {s}(x)} \,\text {d}x + \int _{x_c}^{\Delta x}\! \frac{1}{A_\text {s}(x)} \,\text {d}x \right] ^{-1}. \end{aligned}$$

(34)

For the fluid, the relevant areas are shown in Fig. 14F. The cross-sectional area $A_\text {f}$ is found by subtraction of the solid area from the square that forms a periodic cell

$$\begin{aligned} A_\text {f}(x) = A_{\Box } - A_\text {s}(x), \quad A_{\Box } = 4(\Delta x)^2. \end{aligned}$$

(35)

Again, the corresponding effective thermal transmissibility is found by integration. The integrals were evaluated numerically. Figure 4 shows the computed thermal transmissibility factors (normalised by $\lambda 4 (\Delta x)$) for different fractions of spherical overlap. The normalised thermal transmissibility factor is dimensionless and corresponds to the ratio of effective thermal conductivity and material thermal conductivity, $\lambda _{\mathbb {A}}/ \lambda _\text {s}$, $\lambda _{\mathbb {T}}/ \lambda _\text {f}$, respectively.

Numerical estimation by 3D simulation Since the analytical approach based on the one-dimensional heat equation is not exact, we also perform three-dimensional simulations of a cubic unit cell (with edge length $\Delta x$ due to symmetry, i.e. the unit cell contains 1/8 of a spherical particle). The unit cell is well-resolved with a structured cube grid of ${1\times 10^6}$ cells. The unit cell problem (steady-state heat equation with Dirichlet boundary conditions on left and right boundaries and no-flow otherwise) is solved with a standard finite volume scheme. The results are shown in Fig. 4. We compute effective transmissibility of solid and fluid by evaluating the heat flux over the contact area and the throat, respectively. Effective transmissibilities of the interface are computed by evaluating the heat flux over the void-solid interface. We also compute effective thermal transmissibilities for the entire porous medium by computing flux over the right unit cell face.

Conductive transmissibility model Here we compare results for transmissibilities of the individual phases and the interface from 3D simulation on the cubic unit cell and the model proposed in Sect. 5.2. The model parameters have been fitted manually to provide a good agreement over a wide range of conductivity ratios and particle overlap. The resulting transmissibilities and estimated model parameters are given in Fig. 15.

Fig. 15

Effective thermal conductivities for cubic lattice sintered sphere packing: 3D simulation and model estimation. Normalised solid, void and interface thermal transmissibilities for different fraction of overlap ($\Delta x / R = 0.95$ means 5% overlap with respect to the particle radius) and conductivity ratios. Comparison of 3D unit cell simulations and model fit using the model proposed in Sect. 5.2 with the parameters $C_{\infty ,\text {f}} = A_{\mathbb {T}}A_\text {f}^{-1}$, $C_{0,\text {f}} = 0.1$, $C_{\infty ,\text {s}} = 1.45 A_{\mathbb {A}}A_\text {s}^{-1}$, $C_{0,\text {f}} = 0.45$, $C_{\mathbb {I}}= 0.9$. Areas and distances for the transmissibility model have been computed exactly. The effective areas $A_\text {f}$ and $A_\text {s}$ have been computed as $V_\text {f}/(2\Delta x)$, $V_\text {s}/(2\Delta x)$, respectively, as proposed in Sect. 5.2. The basis for computing transmissibilities from the 3D simulation are the heat fluxes over the fluid throat, the solid contact surface, and the void-solid interface, in combination with the temperature gradients from pore to throat center, grain to contact center, pore to grain center, respectively

Appendix 2

1.1 Convective energy transport at the void-solid interface in wavy channel

Fig. 16

Wavy porous medium geometry and boundary conditions. The upper image shows geometry and partitioning of the wavy channel porous medium. The domain is radially symmetric and rotated around the lower axis. The void-solid interface is described by a sine function with periodicity L. A solid pore is obtained by rotating the plane around the axis by ${90}{^{\circ }}$, a fluid pore is obtained by a full revolution. The bottom image shows boundary conditions for the conjugate heat transfer problem. At the fluid inflow, we prescribe a parabolic velocity profile with given mean velocity $v = v_{\text {f},\text {in}}$ and temperature $T = T_{\text {f},\text {in}}$. At the fluid outflow, the fluid pressure is fixed and heat leaves by advection only. The solid boundary is isolated on inlet and outlet. Over the outer boundary, a constant heat flux is applied. Solid grains occupy quarter sections. The domain boundary faces and void-solid interface faces of grain $S_2$ are shown in green for one exemplary grain in the 3D visualisation

We compute steady-state solutions to the conjugate heat transfer problem (Sect. 2) in a rotational-symmetric wavy channel geometry with DuMu^x. Due to symmetry, a 2D half cross-sectional domain is simulated instead of resolving the full 3D geometry. The fluid flow equations have been adjusted to account for additional terms that appear when formulating the conservative form of the Navier-Stokes equations in cylindrical coordinates (Ferziger and Perić 2008). A 3D render of a section of the channel and the prescribed boundary conditions are shown in Fig. 16. To reduce the influence of the boundary conditions at inlet and outlet on the results, we simulate a channel with six pores and only investigate the two pores in the domain middle. Reynolds numbers have been varied by changing the mean inflow velocity v. The thermal conductivity ratio $\kappa = \lambda _\text {f}/ \lambda _\text {s}$ has been varied by changing $\lambda _\text {s}$. The geometry of the channel has been varied by modifying the pore radius $R_p$ (see Fig. 16). The solid grains are assumed to be cut into four equally shaped grains each resulting from a ${90}{^{\circ }}$ rotation around the symmetry axis. To obtain similar dimensions as for the Berea sandstone sample used in the numerical results section, we used $R_{\mathbb {T}}= {1\times 10^{-5}}\hbox {m}$, $R_p = {2.2\times 10^{-5}}\hbox {m}$, $R_\text {pm} = {5\times 10^{-5}}\hbox {m}$, $L = {6.67\times 10^{-5}}\hbox {m}$ (using the notation of Fig. 16). The interface line is sinusoidal. The resulting porosity is $\phi = 0.22$. Volume-averaged temperatures in pores and grain are computed by integration. The throat temperature is computed as the average of the volume-averaged temperatures of the two adjacent pores. The constant fluid parameters are chosen as $\mu = {1\times 10^{-3}}\hbox {Pa s}$, $\rho = {1000}{\hbox {kg m}^{-3}}$, $\lambda _f = {0.679}{\hbox {W K}^{-1}\hbox {m}^{-1}}$, $c_\text {f}= {4.2}{\hbox {J kg}^{-1} \hbox { K}^{-1}}$, resulting in a Prandtl number of $\text {Pr}= 6.19$.

Fluid flow streamlines for different Reynolds numbers and $\kappa = 0.003$ are shown in Fig. 6. Heat flow lines for two specific cases are shown in Fig. 5. We computed effective thermal conductivities, Nusselt numbers and solid conductivity factors for several geometries, conductivity ratios and Reynolds numbers. The quantities are computed as follows. Let $x_\text {s}$ be the centroid of solid grain S and $x_{\mathbb {T}}$ the centroid of throat ${\mathbb {T}}$ located at the same axial position along the channel. The throat-grain distance is $d_{\mathbb {T}}= \vert \vert x_\text {s}- x_{\mathbb {T}}\vert \vert$. Let $T_{\mathbb {T}}$ denote the average throat temperature and $T_\text {s}$ the average grain temperature (as defined above). We numerically computed $T_{\mathbb {T}}$, $T_\text {s}$, and the mean interface temperature $T_{\mathbb {I}}$ for each solid grain. Local Nusselt numbers and solid conductivity factors are computed by Eqs. (26) and (27),

$$\begin{aligned} \text {Nu}_{\mathbb {T}}= \frac{\vert Q_{{\mathbb {I}}}\vert }{\vert I_{\mathbb {T}}\vert }\frac{d_{\mathbb {T}}}{ \lambda _\text {f}(T_{\mathbb {I}}- T_{\mathbb {T}})}, \quad b_{s,{\mathbb {T}}} = \frac{\vert Q_{{\mathbb {I}}}\vert }{\vert I_{\mathbb {T}}\vert }\frac{d_{\mathbb {T}}}{ \lambda _\text {s}(T_{\mathbb {I}}- T_{K_\text {s}})}, \end{aligned}$$

(36)

where $Q_{{\mathbb {I}}}$ is the total heat flux over the void-solid interface of solid grain $K_\text {s}$. Furthermore, we compute the effective conductivity for convective heat transfer $\lambda _{\text {conv}}$. The results for different conductivity ratios and Reynolds number are shown in Fig. 17. The results for different pore radii and Reynolds are given in Fig. 18. Finally, we suggest a simple empirical mode for estimating the effective thermal conductivity in terms of the Reynolds number. The relation is plotted in Figs. 17 and 18 alongside the simulation results.

Fig. 17

Effective conductivity, Nusselt numbers and solid conductivity factors in wavy channel for different $\kappa$. We considered a wavy channel with six fluid pores. The results are shown for the third and fourth throat in flow direction which is the furthest away from the domain boundary. For small $\kappa < 1$, the heat transfer is best characterised by the Nusselt number and the conductivity factor only plays a minor role. The black solid line shows a simplified heat conductivity model based solely on $\text {Re}_{\mathbb {T}}$. In all cases $\text {Pr}= 6.2$, $R_{\mathbb {T}}= {1\times 10^{-5}}\hbox {m}$, $R_p = {2.2\times 10^{-5}}\hbox {m}$

Fig. 18

Effective conductivity, Nusselt numbers and solid conductivity factors in wavy channel for different pore radii $R_p$. We considered a wavy channel with six fluid pores. The results are shown for the third throat in flow direction which is the furthest away from the domain boundary. Differences in the effective conductivities in the range $1< \text {Re}_{\mathbb {T}}< 50$ can be explained by differences in the local three-dimensional temperature profile. The black solid line shows a simplified heat conductivity model based solely on $\text {Re}_{\mathbb {T}}$. In all cases $\text {Pr}= 6.2$, $\kappa = 0.003$, $R_{\mathbb {T}}= {1\times 10^{-5}}\hbox {m}$

Appendix 3

1.1 Boundary conditions and grid for the micro-cooling scenario

The following section describes in detail how the boundary conditions for the 3D pore-scale conjugate heat transfer simulation and the dual network simulation were chosen to achieve a comparable model scenario.

3D pore-scale conjugate heat transfer simulation We used the ChtMultiRegionFoam module (https://openfoamwiki.net/index.php/ChtMultiRegionFoam) included in OpenFOAM (a conjugate heat transfer solver for the 3D pore-scale equations of Sect. 2) to create numerical reference solutions. We created two Cartesian structured grids with hexagonal axis-aligned cells for the fluid and the solid phase, based on the CT-image voxel data, such that each voxel corresponds to a grid cell. The grid is padded with a fluid layer of one cell thickness $\Delta h = {5.345\times 10^{-6}}\hbox {m}$ on inlet and outlet to ease prescribing boundary conditions. Effectively, this results in fluid and solid having virtually the same temperature at the inlet (${300}{\hbox { K}}$) and the outlet, where a zero temperature gradient condition is prescribed such that $\nabla T_f \cdot {\mathbf {n}}_z = \nabla T_s \cdot {\mathbf {n}}_z = 0$. The bottom plate is realised by padding the grid with a solid layer of one cell thickness, and we assigned the same parameters as for the rest of the solid. The bottom boundary temperature is fixed at ${400}{\hbox { K}}$. Again, this effectively results in fluid and solid having virtually equal temperatures close to the bottom plate. All other domain boundaries are modeled as insulating walls with ${\varvec{v}}_\text {f}= 0$ and $\nabla T_\text {f}\cdot {\mathbf {n}} = \nabla T_\text {s}\cdot {\mathbf {n}} = 0$. Fluid is flowing in positive z-direction by applying fixed pressures at inlet and outlet, proportional to a given pressure gradient.

Dual network simulation The same extracted dual network as described in Sect. 6.1.2 is used for the novel coupled model. For sake of simplicity, we did not add any extra pores or “padding” to the network domain boundaries. Instead, we directly prescribed inlet and outlet pressures (depending on the given pressure gradient) at the fluid pores identified as boundary pores (inlet, outlet) in PoreSpy ’s network extraction algorithm. For the solid and void pores at the bottom plate (same properties as the solid), we directly prescribed $T_\text {f}= T_\text {s}= {400}\hbox { K}$. This can be justified by $\lambda _\text {s}> \lambda _\text {f}$ and follows the observations made in the 3D simulation at the bottom plate boundary where $T_\text {f}$ is very close to ${400}{\hbox { K}}$. The situation is different at the inlet boundary, where a thin padding zone for the fluid is located upstream of the porous domain in the 3D simulation (see Fig. 10). As described above, a fixed inlet temperature of $T_{\mathrm {in}} = {300}{\hbox { K}}$ is applied at the outer boundary of the padding zone in the 3D simulation. For the solid pores at the inlet, we implemented a Robin-type boundary condition for conductive heat flux, taking into account the weighted harmonic average of $\lambda _\text {s}$ and $\lambda _\text {f}$,

$$\begin{aligned} Q_{K_\text {s},\text {in}} = \pi R_p^2 \frac{R_p + \Delta h}{R_p \lambda _\text {s}^{-1} + \Delta h \lambda _\text {f}^{-1}} \frac{T_{\mathrm {in}} - T_{K_\text {s}}}{R_p + \Delta h}, \end{aligned}$$

(37)

where $R_p^2$ is the solid pore’s extended radius and $\Delta h$ the thickness of the 3D simulation padding zone. A similar approach is chosen for the void pores at the inlet:

$$\begin{aligned} Q_{K_\text {f},\text {in}} = \pi R_p^2 \lambda _\text {f}\frac{T_{\mathrm {in}} - T_{K_\text {f}}}{R_p + \Delta h} + \rho c_\text {f}T_{\mathrm {in}} \sum \limits _{{\mathbb {T}}\in {\mathcal {T}}} g_{\mathbb {T}}\left( p_{K_\text {f}} - p_{L_\text {f}}\right) ~. \end{aligned}$$

(38)

Here $R_p$ is the void pore’s inscribed radius. Equation (38) combines heat conduction through the fluid padding zone and convective energy transfer bound to the incoming stream of mass. At the outlet, an outflow condition for the fluid pores was realised by assigning a heat flux

$$\begin{aligned} Q_{K_\text {f},\text {out}} = -\rho c_\text {f}T_{K_\text {f}} \sum \limits _{{\mathbb {T}}\in {\mathcal {T}}} g_{\mathbb {T}}\left( p_{K_\text {f}} - p_{L_\text {f}}\right) \end{aligned}$$

(39)

to leave the domain from control volume $K_\text {f}$. Since the outflow condition implies $\nabla T_\text {s}\cdot {\mathbf {n}} = \nabla T_\text {f}\cdot {\mathbf {n}} = 0$, no outward-bound heat conduction occurs at the solid and void pores located at the outlet boundary. No mass or energy leaves or enters the pores on all other domain boundaries.

Appendix 4

1.1 Energy balance for a micro-cooling device

The main energy fluxes are sketched in Fig. 11. First, we compute the advective energy transport over the inlet face $F_I$ and the outlet face $F_O$

$$\begin{aligned} Q^h_{\text {in,ad}} = \int _{F_I} \rho c_\text {f}T_\text {f}{\varvec{v}}_\text {f}\cdot {\varvec{n}} \text {d}A, \quad Q^h_{\text {out,ad}} = \int _{F_O} \rho c_\text {f}T_\text {f}{\varvec{v}}_\text {f}\cdot {\varvec{n}} \text {d}A. \end{aligned}$$

(40)

The advective inlet heat flux $Q_{h,\text {in}}$ is negative since energy enters the domain by advection, the (advective) outlet heat flux is positive. The inlet fluxes in the pore network are computed by the sum of fluxes into boundary pores. For a single boundary pore, the flux is obtained by

$$\begin{aligned} Q^h_{K_\text {f},\text {in,ad}} = \rho c_\text {f}T_\text {in} \sum \limits _{{\mathbb {T}}\in {\mathcal {T}}} g_{\mathbb {T}}\left( p_{K_\text {f}} - p_{L_\text {f}}\right) , \end{aligned}$$

(41)

where $T_\text {in} = {300}{\hbox { K}}$ and the outlet heat flux is given by Eq. (39). Furthermore, energy enters the system by conduction on the bottom boundary $Q^h_{\text {in,cond}}$ (negative) and leaves by conduction on the inlet $Q^h_{\text {out,cond}}$ (positive). For the 3D simulation, the bottom boundary is solid (padding cells) and the inlet is fluid (padding cells). Hence, the following balances have to be fulfilled

$$\begin{aligned} Q^h_{\text {in,ad}} + Q^h_{\text {out,ad}} + Q^h_{\text {out,cond}} + Q^h_{\text {in,cond}} = 0, \quad Q^h_{\text {in,cond}} + Q^h_{{\mathbb {I}}} = 0, \end{aligned}$$

(42)

with the heat transfer flux

$$\begin{aligned} Q^h_{{\mathbb {I}}} = -\int _{I} \lambda _\text {s}\nabla T_\text {s}\cdot {\varvec{n}} \text {d}A, \end{aligned}$$

(43)

where here I denotes the entire void-solid interface in the sample and ${\varvec{n}}$ is an outer unit normal on the solid surface. For the dual network, since no fluid or solid padding cells are added, the void-solid interface flux does not contain the contributions from the solid boundary at the inlet (outflow of energy by conduction) and the fluid at the bottom boundary (inflow of energy by conduction). These contributions are separately balanced and added to the interface flux to arrive at the same $Q^h_{{\mathbb {I}}}$.