Hostname: page-component-7b7567c8c-jswmv Total loading time: 0 Render date: 2026-11-15T13:10:27.460Z Has data issue: false hasContentIssue false
  • English
  • Français

The first open channel for yield-stress fluids in porous media

Published online by Cambridge University Press:  03 February 2021

Dimitrios FraggedakisOpen the ORCID record for Dimitrios Fraggedakis[Opens in a new window],
Emad ChaparianOpen the ORCID record for Emad Chaparian[Opens in a new window] and
Outi TammisolaOpen the ORCID record for Outi Tammisola[Opens in a new window]
Dimitrios Fraggedakis*
Affiliation:
Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Emad Chaparian
Affiliation:
Mathematics Department, University of British Columbia, Vancouver, BC V6T 1Z2, Canada
Outi Tammisola
Affiliation:
Linné FLOW centre, Department of Engineering Mechanics, KTH Royal Institute of Technology, SE-10044 Stockholm, Sweden
*
†Email address for correspondence: dimfraged@gmail.com
Rights & Permissions [Opens in a new window]

Abstract

The prediction of the first fluidized path of yield-stress fluids in complex porous media is a challenging yet important task to understand the fundamentals of fluid flow in several industrial and biological processes. In most cases, the conditions that open this first path are known either through experiments or expensive computations. Here, we present a simple network model to predict the first open channel for a yield-stress fluid in a porous medium. For porous media made of non-overlapping discs, we find that the pressure drop ${\rm \Delta} P_c$ required to open the first channel for a given yield stress $\tau _y$ depends on both the relative discs size $R_s$ to the macroscopic length $L$ of the system and the packing fraction $\phi$. The non-dimensional pressure gradient ${\rm \Delta} P_c R_s/ \tau _y L$ (i.e. the critical yield number), however, depends on the packing fraction $\phi$ only, leading to a mastercurve for all examined ratios of $R_s/L$. In the case of non-overlapping discs, we find ${\rm \Delta} P_c R_s/ \tau _y L\sim \phi /(1-\phi )$. We also report the statistics on the arclength of the first open path. Finally, we discuss the implication of our results for the design of porous media used in energy storage applications.

JFM classification

Low-Reynolds-number flows: Porous mediaNon-Newtonian Flows: Plastic materials

Information

Type
JFM Papers
Information
Journal of Fluid Mechanics , Volume 911 , 25 March 2021, A58
Check for updates
Copyright
© The Author(s), 2021. Published by Cambridge University Press

1. Introduction

Fluid flows in porous media have been studied for more than a century due to their high relevance to several engineering applications such as enhanced oil recovery (Green & Willhite Reference Green and Willhite1998; Sahimi Reference Sahimi2011; Farajzadeh et al. Reference Farajzadeh, Andrianov, Krastev, Hirasaki and Rossen2012; Fraggedakis et al. Reference Fraggedakis, Kouris, Dimakopoulos and Tsamopoulos2015), filtration and separation (Herzig, Leclerc & Goff Reference Herzig, Leclerc and Goff1970; Tien & Payatakes Reference Tien and Payatakes1979; Jaisi et al. Reference Jaisi, Saleh, Blake and Elimelech2008), fermentation (Pandey Reference Pandey2003; Aufrecht et al. Reference Aufrecht, Fowlkes, Bible, Morrell-Falvey, Doktycz and Retterer2019), soil sequestration (Schlesinger Reference Schlesinger1999), energy storage (Duduta et al. Reference Duduta, Ho, Wood, Limthongkul, Brunini, Carter and Chiang2011; Sun et al. Reference Sun, Zhu, Baumann, Peng, Xu, Shakir, Huang and Duan2019) and food processing (Greenkorn Reference Greenkorn1983). In most cases, the fluids involved in these applications exhibit yield stress/viscoplastic behaviour (Bonn & Denn Reference Bonn and Denn2009; Balmforth, Frigaard & Ovarlez Reference Balmforth, Frigaard and Ovarlez2014). Thus, understanding the conditions – critical pressure drop and/or stresses – that lead to fluidization of yield-stress fluids in porous media can help boost the efficiency and lower the operational cost of several industrial applications.

In pressure-driven flows, the critical pressure drop ${\rm \Delta} P_c$ required to fluidize the yield-stress fluid and open the first channel (Chen, Rossen & Yortsos Reference Chen, Rossen and Yortsos2005; Hewitt et al. Reference Hewitt, Daneshi, Balmforth and Martinez2016; Waisbord et al. Reference Waisbord, Stoop, Walkama, Dunkel and Guasto2019) depends on the heterogeneous geometric characteristics and porosity $1-\phi$ of the porous medium, where $\phi$ is the volume fraction of the solid phase (Talon et al. Reference Talon, Auradou, Pessel and Hansen2013; Bauer et al. Reference Bauer, Talon, Peysson, Ly, Batôt, Chevalier and Fleury2019; Chaparian & Tammisola Reference Chaparian and Tammisola2021). Therefore, it is crucial to understand the relation between the yielding conditions and the structure of the porous medium, which will lead to predictive models for both the first open channel and ${\rm \Delta} P_c$. The classic studies in this era date back to the 1950s and 1960s; the main concern was filtration and oil recovery and limiting pressure gradient models have been proposed (e.g. see Entov Reference Entov1967; Barenblatt, Entov & Ryzhik Reference Barenblatt, Entov and Ryzhik1989). A brief review on the history of the field can be found in Frigaard, Paso & de Souza Mendes (Reference Frigaard, Paso and de Souza Mendes2017).

The classical way to study yield-stress fluids is by solving the momentum equations using viscoplastic constitutive relations, such as the Bingham and Herschel–Bulkley models (Huilgol Reference Huilgol2015; Saramito Reference Saramito2016). The yield-stress behaviour, however, leads to an ill-defined problem that does not describe the stress distribution within the unyielded regions of the fluid (Balmforth et al. Reference Balmforth, Frigaard and Ovarlez2014; Saramito & Wachs Reference Saramito and Wachs2017). Common ways to resolve this problem are by using either optimization- (Hestenes Reference Hestenes1969; Powell Reference Powell1978; Glowinski Reference Glowinski2008; Glowinski & Wachs Reference Glowinski and Wachs2011) or regularization-based methods (Papanastasiou Reference Papanastasiou1987; Frigaard & Nouar Reference Frigaard and Nouar2005). The former are accurate in predicting the yielded/unyielded boundaries and the flow field, however, they are computationally expensive, especially near the yield limit (Saramito & Wachs Reference Saramito and Wachs2017). Although the latter reduce the computational cost, they introduce non-physical parameters that lead to non-physical solutions, incorrect location of yield/unyield boundaries and inaccurate yield limits (Frigaard & Nouar Reference Frigaard and Nouar2005; Tsamopoulos et al. Reference Tsamopoulos, Dimakopoulos, Chatzidai, Karapetsas and Pavlidis2008; Dimakopoulos, Pavlidis & Tsamopoulos Reference Dimakopoulos, Pavlidis and Tsamopoulos2013). Here, we are interested in determining the statistics of the critical ${\rm \Delta} P_c$ for a yield-stress fluid in a porous medium. Thus, we need to use models that can predict accurately and efficiently ${\rm \Delta} P_c$ along with the first open channel.

Fluid flow in a porous medium is traditionally described through network models (Fatt Reference Fatt1956) that represent the complex geometric characteristics of the domain with spherical pore throats and cylindrical edges (Bryant, King & Mellor Reference Bryant, King and Mellor1993; Blunt Reference Blunt2001; Blunt et al. Reference Blunt, Bijeljic, Dong, Gharbi, Iglauer, Mostaghimi, Paluszny and Pentland2013; Alim et al. Reference Alim, Parsa, Weitz and Brenner2017; Stoop et al. Reference Stoop, Waisbord, Kantsler, Heinonen, Guasto and Dunkel2019). In addition to their wide applicability in Newtonian fluids, network models have also been applied to describe ${\rm \Delta} P_c$ and the flow behaviour with respect to the applied pressure drop in yield-stress fluids (Balhoff & Thompson Reference Balhoff and Thompson2004; Balhoff Reference Balhoff2005; Chen et al. Reference Chen, Rossen and Yortsos2005; Sochi Reference Sochi2005; Balhoff et al. Reference Balhoff, Sanchez-Rivera, Kwok, Mehmani and Prodanović2012; Liu et al. Reference Liu, De Luca, Rosso and Talon2019; Talon & Hansen Reference Talon and Hansen2020). When the relation between the local flow rate and the pressure drop is known, the network representation allows for the use of graph theoretic tools (Kharabaf & Yortsos Reference Kharabaf and Yortsos1997; Chen et al. Reference Chen, Rossen and Yortsos2005; Balhoff et al. Reference Balhoff, Sanchez-Rivera, Kwok, Mehmani and Prodanović2012; Liu et al. Reference Liu, De Luca, Rosso and Talon2019) to quickly evaluate ${\rm \Delta} P_c$ and the flow response of the system. In general, however, the results of network viscoplastic models in complex porous media have been rarely compared and validated against those produced by solving the full fluid problem, and thus the conditions of their validity/applicability are still unknown.

The goal of the present work is to predict the first open channel for a yield-stress fluid in a porous medium along with the critical applied pressure drop required to open it. We develop a simple network model based on realistic porous medium configurations, and use graph theoretical tools to study the statistics of the yielding conditions in terms of the medium porosity. We validate our results against pressure-driven simulations of Bingham fluids in porous media. Indeed, we bridge the gap between the continuum flow approaches and network models and also report/analyse the functionality of the critical pressure drop on the porosity. Finally, we discuss the relevance of our study to applications such as semi-solid flow batteries and propose possible extensions.

2. Theory

2.1. Network generation and topology

We are interested in the construction of realistic network models that capture the complex morphology of real porous media. The main scope of our work is to understand the statistics of the critical conditions that lead to fluidization in terms of the porosity $1-\phi$ and the topological characteristics of the medium.

To a first approximation, we assume a porous medium that consists of monodisperse non-overlapping discs (two-dimensional problem) of radius $R_s$, as shown in figure 1(a). It is apparent that the structure of the void space depends on the solid volume fraction defined as $\phi =N_s A_s/A_t$, where $N_s$ is the total number of discs, $A_s = {\rm \pi}R_s^2$ is the area of an individual disc and $A_t=L \times L$ the total area of the system.

Figure 1.

(a) Typical configuration of a porous medium of length $L$ and porosity $1-\phi$ that is made of non-overlapping monodisperse discs of radius $R_s$, and its network representation. Across the domain, a macroscopic pressure gradient ${\rm \Delta} P/L$ is applied. (b) Schematic of three network edges with different $h_i$ and $\ell _i (i=1,2,3)$. The network representation includes the local geometric characteristics of the complex porous medium structure. (c) Characteristic pore-size distribution for the network shown in (a) as derived from the network model. (d) Typical shear stress response $\tau$ as a function of the applied shear rate $\dot {\gamma }$ of a viscoplastic fluid with yield stress $\tau _y$.

For the generation of porous media that consist of non-overlapping discs, we implemented the random sequential addition algorithm (Cule & Torquato Reference Cule and Torquato1999; Torquato, Uche & Stillinger Reference Torquato, Uche and Stillinger2006; Zhang & Torquato Reference Zhang and Torquato2013). The procedure described by Torquato et al. (Reference Torquato, Uche and Stillinger2006) allows for a fast generation of randomly packed discs with the desired volume fractions. Due to the constraint of non-overlapping discs, all generated microstructures never exceed $\phi =0.52$ in two dimensions.

Based on the generated porous medium structure, we can create its network representation as shown in the right image of figure 1(a). The network consists of nodes and edges that span the entire medium, where its complex topological characteristics are encoded on the connectivity between them (Gostick Reference Gostick2017). The local geometric characteristics of the porous medium are included in the length $\ell _i$ and gap $2h_i$ of each individual edge, as we show in figure 1(b).

For the generation of the network model, we implemented the maximal ball algorithm (Silin & Patzek Reference Silin and Patzek2006; Al-Kharusi & Blunt Reference Al-Kharusi and Blunt2007; Dong & Blunt Reference Dong and Blunt2009) which allows us to obtain the $h_i$ and length $\ell _i$ for each edge. The maximal ball algorithm ‘fits’ a circle/sphere within each pore (Alim et al. Reference Alim, Parsa, Weitz and Brenner2017), the radius of which represents $h_i$. As will be discussed in § 3.2, there are recently developed methods (Gostick Reference Gostick2017) that can extract $h_i$ and $\ell _i$ using image processing, such as the watershed segmentation. Other choices of $h_i$ can be considered (e.g. equivalent gap/radius etc.), however, our choice for the local radius works well in predicting both the first open channel and the critical pressure drop required to open it. The generated network was represented by a graph with vertices $V$ and edges $E$ using the open source library NetworkX (Hagberg, Swart & Chult Reference Hagberg, Swart and Chult2008). For demonstration, we show in figure 1(c) the pore-size distribution for the configuration of figure 1(a).

2.2. Yield-stress fluid in a network

Yield-stress fluids are characterized by their solid–liquid transition when the Euclidean norm of the stress field exceeds the value of yield stress (i.e. von Mises criterion; see Hill Reference Hill1998; Gurtin, Fried & Anand Reference Gurtin, Fried and Anand2010). The typical shear stress response in simple shear flow is shown in figure 1(d), where for $\dot {\gamma }\rightarrow 0$ the shear stress reaches its critical value $\tau \rightarrow \tau _y$. The most common constitutive relations used to describe ‘simple’ yield-stress fluids are the Bingham and Herschel–Bulkley models, which differ in their viscosity functional. Both of them, however, predict the same yield limit (Frigaard Reference Frigaard2019) since the viscous dissipation is negligible compared to the plastic one for $\dot {\gamma }\rightarrow 0$. Therefore, it is sufficient to discuss only the Bingham model for a porous medium to understand the connection between the yielding conditions and the geometric and topological characteristics of the network.

For a Bingham fluid pressure-driven flow close to the yield limit (i.e. the linearization of the flow rule), the local flow rate $q_i$ of edge $i$ is described in terms of the local geometric properties $h_i,\ell _i$, the local pressure drop along the edge ${\rm \Delta} P_i$ and the yield stress $\tau _y$ of the fluid as,

(2.1)\begin{equation} q_{i} \sim \left\{ \begin{array}{@{}ll} \dfrac{h_i^3}{\ell_i} \left({\rm \Delta} P_{i}-\dfrac{\tau_y \ell_i}{h_i}\right) & \text{for } {\rm \Delta} P_{i}>\dfrac{\tau_y \ell_i}{h_i},\\ 0 & \text{for } \lvert{\rm \Delta} P_{i}\rvert < \dfrac{\tau_y \ell_i}{h_i}, \\ \dfrac{h_i^3}{\ell_i} \left({\rm \Delta} P_{i}+\dfrac{\tau_y \ell_i}{h_i}\right) & \text{for } {\rm \Delta} P_{i}<-\dfrac{\tau_y \ell_i}{h_i}.\\ \end{array} \right. \end{equation}

Near the no-flow limit $q_i\rightarrow 0$, we see from (2.1) that ${\rm \Delta} P_i\rightarrow \tau _y\ell _i/h_i$, and thus channels that are smaller in radius and/or longer in length require a larger applied pressure drop to yield. This is the common yield criterion for plane Poiseuille flow of a Bingham fluid in a channel of width $2h_i$ and length $\ell _i$. Across the first open channel, we can calculate the total pressure drop across the medium to be ${\rm \Delta} P_c \equiv \sum _{i=1}^N {\rm \Delta} P_{c,i} = \tau _y\sum _{i=1}^N \ell _i/h_i$, where $N$ is the total number of edges across the path. From this expression, we can see that the connectivity between the edges determines the first open channel in a real porous medium, and it corresponds to the path of ‘least resistance’. Thus, the problem of finding ${\rm \Delta} P_c$ can be formulated as finding the path of the minimum pressure drop as follows (Liu et al. Reference Liu, De Luca, Rosso and Talon2019)

(2.2)\begin{equation} \frac{{\rm \Delta} P_c}{\tau_y} = \min\limits_{C \in \mathcal{C}_{{in-out}}} \sum_{i=1}^N \frac{\ell_i}{h_i}, \end{equation}

where $\mathcal {C}_{{in-out}}$ is the set of all paths between the corresponding boundaries. The problem of (2.2) satisfies the principle of minimum dissipation rate and is valid near equilibrium (De Groot & Mazur Reference De Groot and Mazur1984). In particular, for isothermal pressure-driven flow of temperature $T$, the energy dissipation according to linear irreversible thermodynamics is $\varPhi =q{\rm \Delta} P$ (De Groot & Mazur Reference De Groot and Mazur1984). Thus, for conditions near the solid–liquid transition where $q\rightarrow 0^+$, the minimum pressure drop path also minimizes $\varPhi$.

To solve equation (2.2), we transform the generated network into a graph with edges that have weights equal to $\ell _i/h_i$ and use the Dijkstra method (Dijkstra Reference Dijkstra1959) for directed graphs to determine the first open channel. This method is known to scale quadratically with the path length (Bollobás Reference Bollobás2013), and therefore the computational cost increases for complex domains with a larger number of edges. For a single porous medium configuration, however, the overall computational time to determine the first open channel is much lower (seconds to minutes) than that required to solve the full fluid flow problem using optimization methods (days to weeks) (Dimakopoulos et al. Reference Dimakopoulos, Makrigiorgos, Georgiou and Tsamopoulos2018; Chaparian & Tammisola Reference Chaparian and Tammisola2021).

3. Results

3.1. The first open channel

When the applied pressure drop approaches the critical value ${\rm \Delta} P \rightarrow {\rm \Delta} P_c^+$, there exists a single open channel across the entire medium. Here, we test the validity of the proposed approach to determine the first open channel when the solid–liquid transition occurs. For comparison, we solve the full flow field under pressure-driven conditions for a yield-stress fluid for the porous media shown in figure 2. We consider the cases of $\phi =0.3$ and $\phi =0.5$, respectively. All the lengths are normalized with the macroscopic length of the system $L$ and also $R_s/L=0.02$ and $R_s/L=0.1$. In this study, we investigate two-dimensional porous media, however, the network approach is general and can also be used to study three-dimensional geometries.

Figure 2.

Validation of the network model for the first open channel against simulations for a pressure-driven Bingham fluid in a porous medium with $R_s/L=0.02$ and $R_s/L=0.1$. (a) Simulation results of the open pathway for conditions near the critical pressure drop ${\rm \Delta} P_c$. The contour plot shows the magnitude of the local velocity, normalized with the maximum velocity across the channel. (b) Network model predictions for the first open channel. The cases of $\phi =0.3$ and $\phi =0.5$ are examined. It is clear that both the full Bingham fluid simulation and the network model predict the same location for the first open channel.

Figure 2(a) shows the normalized velocity magnitude from the solution of the momentum equation (i.e. Stokes equation) for a Bingham fluid (Chaparian & Tammisola Reference Chaparian and Tammisola2021). Details of the numerical simulation of the fluid flow problem shown in figure 2(a) are discussed by Chaparian & Tammisola (Reference Chaparian and Tammisola2021). To highlight, the porous medium is designed with a desirable porosity by randomizing non-overlapping monodisperse discs. To solve the equations of motion (Stokes equation), we use the augmented Lagrangian method implemented in a finite element environment (FreeFEM++; see Hecht Reference Hecht2012) coupled with adaptive mesh (Roquet & Saramito Reference Roquet and Saramito2003). At the inlet (left side of the computational box) we set the flow rate (see Roustaei & Frigaard (Reference Roustaei and Frigaard2015), for the algorithm), and as a result, we calculate the pressure drop for different values of yield stress. More precisely, the velocity is scaled with the average inlet velocity $U$, hence, the non-dimensional flow rate is always equal to unity. Therefore, the yield limit is transferred to $\tau _y \ell _{ch} / \mu U \to \infty$ (i.e. Bingham number $\to \infty$), where $\ell _{ch}$ is a characteristic length. When $\tau _y \ell _{ch} / \mu U \to \infty$ the non-dimensional pressure drop goes to infinity as well ${\rm \Delta} P \ell _{ch}^2 / \mu U L \to \infty$. However, the ratio of the non-dimensional pressure drop to the Bingham number, ${{\rm \Delta} P \ell _{ch}}/{\tau _y L}$ or the yield number, is of interest, which at this limit asymptotically reaches a finite value. So the limiting pressure drop to open the first channel can be computed by considering the limit $\tau _y \ell _{ch} / \mu U \to \infty$. For more details, we refer the interested reader to Chaparian et al. (Reference Chaparian, Izbassarov, De Vita, Brandt and Tammisola2020) and Chaparian & Tammisola (Reference Chaparian and Tammisola2021).

Figure 2(b) depicts the predictions for the first open channel after solving the minimization problem of (2.2). In all cases (i.e. $\phi =0.3$ and $0.5$ and $R_s/L=0.02$ and $0.1$), the network model is able to reproduce the results of the fluid flow problem (i.e. (a)) for the first open path. In the full fluid flow problem for $R_s/L=0.02$, we can see the existence of minor tiny ‘roundabout’ paths in the vicinity of the major path. However, it is clear from figure 2(a) that the flow rate in those minor paths is negligible and does not contribute to the leading order of the dissipation that is attributed to the major route predicted by the network model, figure 2(b).

3.2. Critical pressure drop ${\rm \Delta} P_c$ and its statistics in porous media

In addition to the first open path, the network model can predict the critical pressure difference ${\rm \Delta} P_c$ required to yield the fluid in the porous medium.

In table 1 we show the predictions of the normalized critical pressure drop ${\rm \Delta} P_c/\tau _y$ for both the network model and the full simulations. For $\phi =0.3$ and $0.5$ we consider the configurations shown in figure 2. In the cases of $R_s/L=0.02$, it is clear that the relative difference in the predicted ${\rm \Delta} P_c$ by the two approaches never exceeds $6.5\,\%$. The negligible difference can be justified by the simplification of the geometric characteristics of the medium by its network representation. For $R_s/L=0.1$, the relative error is approximately $12.1\,\%$, which we attribute to the large curvature effects of the discs that the network model cannot capture. We conclude, however, that network-based models are adequate to predict both the first open path and the critical pressure drop required to open it.

Table 1.

Predicted normalized critical pressure drop ${\rm \Delta} P_c$ required to open the first open channel for the configurations shown in figure 2. We show both the predictions of the network model and the full simulation results, along with the relative error $({\rm \Delta} P_{c}^{{sim}}-{\rm \Delta} P_{c}^{{net}} ) / {\rm \Delta} P_{c}^{{sim}}$. The network model is adequate on predicting both the first open channel and the critical pressure drop required to open it.

Before we proceed with the analysis of the model, we want to comment on the choice of $h_i$ and $\ell _i$ for calculating the critical pressure drop. Given the geometric irregularities of porous media, the choice of the local geometric quantities present in the model is non-trivial. Here, in the case of non-overlapping discs, we used the $h_i$ and $\ell _i$ that were directly calculated by the maximum-ball algorithm. Other choices are possible, where for instance one can use the equivalent gap length for evaluating $h_i$ defined as the ratio of the area of the channel to its perimeter. To better understand the influence of such a choice, we repeated the calculations for $\phi =0.5$ and $R_s/L=0.02$ where now $h_i$ is the equivalent gap length and has been extracted by the image processing methods described by Gostick et al. (Reference Gostick, Khan, Tranter, Kok, Agnaou, Sadeghi and Jervis2019). While the results for the first open path are identical, we find ${\rm \Delta} P_{c,equiv}/\tau _y\simeq 1016.6$, which is an order of magnitude larger than the fluid flow simulation predictions. This comparison indicates that the critical pressure drop is sensitive to the choice of $h_i$ and $\ell _i$, however, the location as well as the topology of the first open path tend to be unchanged.

The validation of the model allows us to predict ${\rm \Delta} P_c$ as a function of $\phi$ and the geometric characteristics of the system. For a porous medium made of non-overlapping monodispersed discs, we control the microstructure characteristics by changing the ratio $R_s/L$. For each combination of $\phi$ and $R_s/L$, we generate 500 realizations to gather the statistics of ${\rm \Delta} P_c$.

Figure 3(a) shows the normalized critical pressure drop ${\rm \Delta} P_c/\tau _y$ in terms of $\phi$. The coloured lines indicate different values of $R_s/L\in [0.02,0.1]$, while the error bars correspond to the variance of the statistical sample. For all $R_s/L$ the normalized pressure drop increases with increasing $\phi$. This behaviour is expected as $h_i$ of each edge decreases monotonically as $\phi$ increases (Torquato & Haslach Reference Torquato and Haslach2002). Additionally, we observe that decreasing the ratio $R_s/L$ leads to an increase in ${\rm \Delta} P_c$, in qualitative agreement with experiments (Waisbord et al. Reference Waisbord, Stoop, Walkama, Dunkel and Guasto2019). The reason for this behaviour can be explained by examining the arclength (defined as $L_c=\sum _i^N \ell _i$) for the first open channel.

Figure 3.

Statistics of the predictions of the network model for 500 realizations per value of $\phi$. (a) Critical pressure drop ${\rm \Delta} P_c$ as a function of the volume fraction $\phi$ for different ratios of $R_s/L$; ${\rm \Delta} P_c$ is normalized with the yield stress of the fluid $\tau _y$. The error bars indicate the variance around the mean value. Inset – scaled critical pressure drop ${\rm \Delta} P_c R_s/\tau _y L$ in terms of $\phi$. For non-overlapping discs, the results for different $R_s/L$ collapse into a master curve. (b) Probability density distributions for the normalized arclength of the first open channel $L_c/L-1$. The cases of $R_s/L=0.02$ and $R_s/L=0.08$ are shown, respectively, for $\phi =0.1$ and $\phi =0.5$. In general, increasing $\phi$ leads to increase of ${\rm \Delta} P_c$ required to open the first channel, but also leads to a wide distribution of arclengths, which provide a large uncertainty in the critical macroscopic pressure gradient ${\rm \Delta} P_c/L_c$ required to yield the fluid in the porous medium.

Figure 3(b) illustrates the histogram of $L_c$ for $\phi =0.3$ and $0.5$ as well as $R_s/L=0.02$ and $0.08$. In both cases we observe that increasing the solid volume fraction leads to an increase in the total length of the first open path. While large $\phi$ results in almost similar behaviour for both $R_s/L$, it is apparent that decreasing the size of the discs results in more tortuous paths, even for low solid volume fraction.

Dimensional analysis of (2.2) indicates that ${\rm \Delta} P_c$ follows a simple scaling with $R_s/L$. In particular, by rescaling the local $h_i$ with the discs radius $R_s$, we find ${\rm \Delta} \tilde {P}_c={\rm \Delta} P_c R_s/\tau _y L = \sum _{i=1}^N (\ell _i/L)(R_s/h_i)$. This form can be interpreted as the ratio between the applied pressure drop ${\rm \Delta} P/L$ over the local shear force per volume generated by the existence of particles with radius $R_s$.

The inset in figure 3(a) shows the rescaled form of the critical pressure drop for all $R_s/L$, where for non-overlapping discs a mastercurve exists for all the examined values of $\phi$. To better understand the functional dependence of ${\rm \Delta} \tilde {P}_c$ with $\phi$, we follow a mean-field approach and approximate its expression as

(3.1)\begin{equation} {\rm \Delta} \tilde{P}_c\simeq N \left\langle \frac{\ell}{L} \right\rangle \left\langle \frac{R_s}{h} \right\rangle = \left\langle \frac{L_c}{L} \right\rangle \left\langle \frac{R_s}{h} \right\rangle, \end{equation}

where $N$ is the number of segments of the path, while $\langle \ell \rangle$ and $\langle 1/h \rangle$ are the mean length and inverse gap of each edge, respectively. To find the functional dependence for both $\langle {L_c}/{L} \rangle$ and $\langle {R_s}/{h} \rangle$, we calculate their average directly from our collected data values where we find $\langle {L_c}/{L} \rangle \sim 0.5\phi$ and $\langle {R_s}/{h} \rangle \sim 2/(1-\phi )$; see figures 4(a) and 4(b). In the porous medium community, $L_c/L$ is known as the geometric tortuosity $\tau _g$ and several analytical and empirical relations have been proposed to describe its form as a function of the porosity $(1-\phi )$ (Guo Reference Guo2012). To a first approximation, we describe $\langle L_c/L \rangle$ with the Bruggeman relation for cylinders $\langle L_c/L \rangle \sim (1-\phi )^{-1/2}$ (Tjaden et al. Reference Tjaden, Cooper, Brett, Kramer and Shearing2016). In the limit of $\phi \rightarrow 0$, we find $\langle L_c/L \rangle \simeq 1+0.5\phi + O(\phi ^2)$, which shows a very similar behaviour to the scaling reported in figure 4(b). While this limit corresponds to a dilute mixture of particles in a permeable matrix, the geometric characteristics of the first open roughly follow this estimation in the range of studied porosities.

Figure 4.

Scalings of the geometric characteristics of a porous medium made of non-overlapping discs. (a) Inverse channel gap $\langle R_s/h \rangle$ thickness vs. $(1-\phi )^{-1}$. From our calculations it is clear that $\langle R_s/h \rangle$ scales with the inverse of the porosity $1-\phi$. (b) First open channel length $\langle L_c/L \rangle$ vs. $\phi$. Based on a simple estimate of the geometric tortuosity, we find that $\langle L_c/L \rangle$ scales linearly with $\phi$. (c) Non-dimensional pressure gradient ${\rm \Delta} \tilde {P}_c={\rm \Delta} P_c R_s/\tau _y L$ vs. $\phi /(1-\phi )$. According to (3.1), the product of $\langle L_c/L \rangle$ and $\langle R_s/h \rangle$ indicates that ${\rm \Delta} \tilde {P}_c$ scales with $\phi /(1-\phi )$.

To this end, we can approximate ${\rm \Delta} \tilde {P}_c\simeq Af(\phi )$ for monodisperse discs of radius $R_s$, where $A$ is a proportionality constant and $f(\phi )=\phi /(1-\phi )$. In figure 4(c) we show the dimensionless pressure drop ${\rm \Delta} \tilde {P}_c$ vs. $f(\phi )$, where we see a linear relationship between the two quantities, with $A\simeq 3$, validating the mean-field assumption considered in (3.1).

4. Summary and discussion

4.1. Summary

In this work, we investigated the application of network models for yield-stress fluids in porous media. We find that network models provide an efficient and accurate way to describe the fluidization conditions in porous media. In particular, they capture the first open channel, which is equivalent to the path of least resistance through the entire medium. We also demonstrated the capabilities of the network model to predict the critical applied pressure drop required to open the first channel in the medium. This was rationalized by comparing the network results with direct numerical simulations of the fluid flow problem. We showed the accuracy and computational efficiency of the network-based models in the context of solid–liquid transition. Also, we performed statistical analyses on different two-dimensional porous media made of non-overlapping discs to understand the dependence of the critical pressure drop ${\rm \Delta} P_c$ on the solid volume fraction $\phi$ and the geometric characteristics of the solid matrix. For two-dimensional porous media that consist of non-overlapping discs, we revealed that the ratio ${\rm \Delta} P_c R_s/\tau _y L$ scales linearly with $\phi /(1-\phi )$, which is a result of the average geometric characteristics of the first open path.

4.2. Discussion

The present model considers only the case of viscoplastic materials and disregards the fluid elasticity (Saramito Reference Saramito2007; Fraggedakis, Dimakopoulos & Tsamopoulos Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a,Reference Fraggedakis, Dimakopoulos and Tsamopoulosb; Chaparian et al. Reference Chaparian, Izbassarov, De Vita, Brandt and Tammisola2020) prior to yielding. Therefore, further analysis is required to connect the yield criterion to the elastic modulus of the fluid, which can provide insights for engineering both the porous medium and the fluid itself. Additionally, the yielding and/or stoppage conditions might be further affected by possible thixotropic (Mewis & Wagner Reference Mewis and Wagner2012) and kinematic hardening (Gurtin et al. Reference Gurtin, Fried and Anand2010; Dimitriou & McKinley Reference Dimitriou and McKinley2019) phenomena.

Beyond the rheological characteristics of the model, more complete closures for the yielding criterion, (2.2), can be considered that take into account the geometric irregularities of the channels. For instance, Roustaei et al. (Reference Roustaei, Chevalier, Talon and Frigaard2016) showed that ${\rm \Delta} P_{c,i}$, even for quite smooth channels, can be up to 25 % different compared to that predicted for a typical planar Poiseuille flow. Also, irregular channels that show fracture-like characteristics tend to yield at a smaller applied pressure drop compared to straight channels. Another complication on the yield-stress criterion might come from the dimensionality of the porous medium. In three-dimensional structures, the description of open channels as tubes might not be adequate, and the geometric quantities $h_i$ and $\ell _i$ need to be chosen carefully (Waisbord et al. Reference Waisbord, Stoop, Walkama, Dunkel and Guasto2019). Thus, for irregular porous media, the yielding criterion, and thus the expression related to the minimum pressure to open the first channel, should be modified.

The mean-field approximation for the critical pressure drop ${\rm \Delta} \tilde {P}_c$ as given in (3.1), indicates that one can develop analytical approximations for the critical pressure drop based only on geometries quantities of a medium that consists of non-overlapping objects. In particular, one can follow a more rigorous approach to estimate $\langle R_s/h \rangle$ analytically, as presented by Torquato & Haslach (Reference Torquato and Haslach2002). In particular, $\langle R_s/h \rangle =\langle 1/(l-1) \rangle$ where $l$ is the nearest-neighbour distance between particles, and can be formally calculated as

(4.1)\begin{equation} \left\langle \frac{1}{l-1} \right\rangle = \int_{1}^{\infty} \frac{1}{l-1}H_P(l)\, \textrm{d} l, \end{equation}

where $H_P$ is the ‘particle’ nearest-neighbour probability density function and depends on the geometry and dimensionality of the particles (Torquato, Lu & Rubinstein Reference Torquato, Lu and Rubinstein1990) as well as on the protocol used to generate the particle configuration. The random sequential addition algorithm used here produces non-equilibrium configurations where no analytical expressions for $H_P$ are known. However, when equilibrium configurations are generated (e.g. through Monte Carlo algorithms), Torquato (Reference Torquato1995) has reported exact and approximate expressions for $H_P$ for monodisperse rods, discs and spheres.

Our results on the normalized critical pressure drop can be used to design porous media systems with the desired flow properties. In applications such as semi-solid flow batteries (Duduta et al. Reference Duduta, Ho, Wood, Limthongkul, Brunini, Carter and Chiang2011), it is critical to keep the contact between the active material (electrode particles) and the conductive wiring (carbon nanoparticle network) intact during operation (Wei et al. Reference Wei, Fan, Helal, Smith, McKinley, Chiang and Lewis2015). This can be achieved by immersing the active material and the electronically conductive agent in yield-stress fluids like Carbopol (Zhu et al. Reference Zhu2020). Therefore, we can use the predicted ${\rm \Delta} P_c/\tau _y$ to determine the size of the active particles to optimize the design and operation of semi-solid electrodes.

The present model can also provide insights into the design of porous media. By taking advantage of the computational efficiency of the proposed network model, we can perform on-the-fly optimization to construct porous media with optimal mixing and transport properties (Lester, Metcalfe & Trefry Reference Lester, Metcalfe and Trefry2013). Such ideas have recently been implemented in elastic networks with optimal phonon band structures (Ronellenfitsch & Dunkel Reference Ronellenfitsch and Dunkel2019), and we believe they can also be used for designing porous media immersed in yield-stress fluids.

Given their inherent node/edge structure, network models are fairly simple to analyse using graph theoretical tools. The unique property of yield-stress fluid, namely the fluidization condition, allows us to use algorithms that can find the minimum resistance pathways with minimal effort. In coarse-grained domains, however, where the microscopic geometric irregularities are encoded in the heterogeneous ‘permeability’ tensor (Hewitt et al. Reference Hewitt, Daneshi, Balmforth and Martinez2016), graph theory tools might not be the most suitable ones. An alternative way to calculate the first open channel in a continuum with spatially variable properties is through methodologies used in physical chemistry to identify reaction pathways (Henkelman, Uberuaga & Jónsson Reference Henkelman, Uberuaga and Jónsson2000). There, the first open channel corresponds to the path that passes through the minimum energy barrier, namely the transition state point.

Acknowledgements

D.F. (aka dfrag) wants to thank T. Zhou, M. Mirzadeh, A. Sahu and Y. Dimakopoulos for insightful discussions related to the validity of the network model.

Declaration of interest

The authors report no conflict of interest.

Author contributions

D.F. conceptualized and designed the present study. D.F. and E.C. performed the analysis and contributed equally in the present work. D.F. wrote the first version of the present manuscript. All authors contributed to the final manuscript.

References

REFERENCES

Al-Kharusi, A.S. & Blunt, M.J. 2007 Network extraction from sandstone and carbonate pore space images. J. Petrol. Sci. Engng 56 (4), 219–231.CrossRefGoogle Scholar
Alim, K., Parsa, S., Weitz, D.A. & Brenner, M.P. 2017 Local pore size correlations determine flow distributions in porous media. Phys. Rev. Lett. 119 (14), 144501.CrossRefGoogle ScholarPubMed
Aufrecht, J.A., Fowlkes, J.D., Bible, A.N., Morrell-Falvey, J., Doktycz, M.J. & Retterer, S.T. 2019 Pore-scale hydrodynamics influence the spatial evolution of bacterial biofilms in a microfluidic porous network. PLoS ONE 14 (6), e0218316.CrossRefGoogle Scholar
Balhoff, M. 2005 Modeling the flow of non-Newtonian fluids in packed beds at the pore scale. PhD thesis, Louisiana State University.Google Scholar
Balhoff, M., Sanchez-Rivera, D., Kwok, A., Mehmani, Y. & Prodanović, M. 2012 Numerical algorithms for network modeling of yield stress and other non-Newtonian fluids in porous media. Transp. Porous Med. 93 (3), 363–379.CrossRefGoogle Scholar
Balhoff, M.T. & Thompson, K.E. 2004 Modeling the steady flow of yield-stress fluids in packed beds. AIChE J. 50 (12), 3034–3048.CrossRefGoogle Scholar
Balmforth, N.J., Frigaard, I.A. & Ovarlez, G. 2014 Yielding to stress: recent developments in viscoplastic fluid mechanics. Annu. Rev. Fluid Mech. 46, 121–146.CrossRefGoogle Scholar
Barenblatt, G.I., Entov, V.M. & Ryzhik, V.M. 1989 Theory of Fluid Flows through Natural Rocks. Kluwer Academic Publishers.Google Scholar
Bauer, D., Talon, L., Peysson, Y., Ly, H.B., Batôt, G., Chevalier, T. & Fleury, M. 2019 Experimental and numerical determination of Darcy's law for yield stress fluids in porous media. Phys. Rev. Fluids 4 (6), 063301.CrossRefGoogle Scholar
Blunt, M.J. 2001 Flow in porous media—pore-network models and multiphase flow. Curr. Opin. Colloid Interface Sci. 6 (3), 197–207.CrossRefGoogle Scholar
Blunt, M.J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A. & Pentland, C. 2013 Pore-scale imaging and modelling. Adv. Water Resour. 51, 197–216.CrossRefGoogle Scholar
Bollobás, B. 2013 Modern Graph Theory, vol. 184. Springer Science & Business Media.Google Scholar
Bonn, D. & Denn, M.M. 2009 Yield stress fluids slowly yield to analysis. Science 324 (5933), 1401–1402.CrossRefGoogle Scholar
Bryant, S.L., King, P.R. & Mellor, D.W. 1993 Network model evaluation of permeability and spatial correlation in a real random sphere packing. Transp. Porous Med. 11 (1), 53–70.CrossRefGoogle Scholar
Chaparian, E., Izbassarov, D., De Vita, F., Brandt, L. & Tammisola, O. 2020 Yield-stress fluids in porous media: a comparison of viscoplastic and elastoviscoplastic flows. Meccanica 55 (2), 331–342.CrossRefGoogle ScholarPubMed
Chaparian, E. & Tammisola, O. 2021 Sliding flows of yield-stress fluids. Accepted for publication in J. Fluid Mech. (preprint) arXiv:2004.02950.CrossRefGoogle Scholar
Chen, M., Rossen, W. & Yortsos, Y.C. 2005 The flow and displacement in porous media of fluids with yield stress. Chem. Engng Sci. 60 (15), 4183–4202.CrossRefGoogle Scholar
Cule, D. & Torquato, S. 1999 Generating random media from limited microstructural information via stochastic optimization. J. Appl. Phys. 86 (6), 3428–3437.CrossRefGoogle Scholar
De Groot, S.R. & Mazur, P. 1984 Non-Equilibrium Thermodynamics. Dover Publication Inc.Google Scholar
Dijkstra, E.W. 1959 A note on two problems in connexion with graphs. Numer. Math. 1 (1), 269–271.CrossRefGoogle Scholar
Dimakopoulos, Y., Makrigiorgos, G., Georgiou, G.C. & Tsamopoulos, J. 2018 The PAL (Penalized Augmented Lagrangian) method for computing viscoplastic flows: a new fast converging scheme. J. Non-Newtonian Fluid Mech. 256, 23–41.CrossRefGoogle Scholar
Dimakopoulos, Y., Pavlidis, M. & Tsamopoulos, J. 2013 Steady bubble rise in Herschel–Bulkley fluids and comparison of predictions via the augmented Lagrangian method with those via the Papanastasiou model. J. Non-Newtonian Fluid Mech. 200, 34–51.CrossRefGoogle Scholar
Dimitriou, C.J. & McKinley, G.H. 2019 A canonical framework for modeling elasto-viscoplasticity in complex fluids. J. Non-Newtonian Fluid Mech. 265, 116–132.CrossRefGoogle Scholar
Dong, H. & Blunt, M.J. 2009 Pore-network extraction from micro-computerized-tomography images. Phys. Rev. E 80 (3), 036307.CrossRefGoogle ScholarPubMed
Duduta, M., Ho, B., Wood, V.C., Limthongkul, P., Brunini, V.E., Carter, W.C. & Chiang, Y. 2011 Semi-solid lithium rechargeable flow battery. Adv. Energy Mater. 1 (4), 511–516.CrossRefGoogle Scholar
Entov, V.M. 1967 On some two-dimensional problems of the theory of filtration with a limiting gradient. Prikl. Mat. Mekh. 31 (5), 820–833.Google Scholar
Farajzadeh, R., Andrianov, A., Krastev, R., Hirasaki, G.J. & Rossen, W.R. 2012 Foam–oil interaction in porous media: implications for foam assisted enhanced oil recovery. Adv. Colloid Interface Sci. 183, 1–13.CrossRefGoogle ScholarPubMed
Fatt, I. 1956 The network model of porous media. Trans. AIME 207 (1), 144–181.CrossRefGoogle Scholar
Fraggedakis, D., Dimakopoulos, Y. & Tsamopoulos, J. 2016 a Yielding the yield-stress analysis: a study focused on the effects of elasticity on the settling of a single spherical particle in simple yield-stress fluids. Soft Matt. 12 (24), 5378–5401.CrossRefGoogle ScholarPubMed
Fraggedakis, D., Dimakopoulos, Y. & Tsamopoulos, J. 2016 b Yielding the yield stress analysis: a thorough comparison of recently proposed elasto-visco-plastic (EVP) fluid models. J. Non-Newtonian Fluid Mech. 236, 104–122.CrossRefGoogle Scholar
Fraggedakis, D., Kouris, C., Dimakopoulos, Y. & Tsamopoulos, J. 2015 Flow of two immiscible fluids in a periodically constricted tube: transitions to stratified, segmented, churn, spray, or segregated flow. Phys. Fluids 27 (8), 082102.CrossRefGoogle Scholar
Frigaard, I.A. 2019 Background lectures on ideal visco-plastic fluid flows. In Lectures on Visco-Plastic Fluid Mechanics, pp. 1–40. Springer.CrossRefGoogle Scholar
Frigaard, I.A. & Nouar, C. 2005 On the usage of viscosity regularisation methods for visco-plastic fluid flow computation. J. Non-Newtonian Fluid Mech. 127 (1), 1–26.CrossRefGoogle Scholar
Frigaard, I.A., Paso, K.G. & de Souza Mendes, P.R. 2017 Bingham's model in the oil and gas industry. Rheol. Acta 56 (3), 259–282.CrossRefGoogle Scholar
Glowinski, R. 2008 Lectures on Numerical Methods for Non-Linear Variational Problems. Springer Science & Business Media.Google Scholar
Glowinski, R. & Wachs, A. 2011 On the numerical simulation of viscoplastic fluid flow. In Handbook of Numerical Analysis, vol. 16, pp. 483–717. Elsevier.CrossRefGoogle Scholar
Gostick, J.T. 2017 Versatile and efficient pore network extraction method using marker-based watershed segmentation. Phys. Rev. E 96 (2), 023307.CrossRefGoogle ScholarPubMed
Gostick, J.T., Khan, Z.A., Tranter, T.G., Kok, M.D.R., Agnaou, M., Sadeghi, M. & Jervis, R. 2019 Porespy: a python toolkit for quantitative analysis of porous media images. J. Open Source Softw. 4 (37), 1296.CrossRefGoogle Scholar
Green, D.W. & Willhite, G.P. 1998 Enhanced Oil Recovery, vol. 6. Society of Petroleum Engineers.Google Scholar
Greenkorn, R.A. 1983 Flow Phenomena in Porous Media: Fundamentals and Applications in Petroleum, Water and Food Production. Marcel Dekker.Google Scholar
Guo, P. 2012 Dependency of tortuosity and permeability of porous media on directional distribution of pore voids. Transp. Porous Med. 95 (2), 285–303.CrossRefGoogle Scholar
Gurtin, M.E., Fried, E. & Anand, L. 2010 The Mechanics and Thermodynamics of Continua. Cambridge University Press.CrossRefGoogle Scholar
Hagberg, A., Swart, P. & Chult, D.S. 2008 Exploring network structure, dynamics, and function using NetworkX. Tech. Rep. Los Alamos National Laboratory (LANL).Google Scholar
Hecht, F. 2012 New development in freefem++. J. Numer. Math. 20 (3), 251–265.CrossRefGoogle Scholar
Henkelman, G., Uberuaga, B.P. & Jónsson, H. 2000 A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 113 (22), 9901–9904.CrossRefGoogle Scholar
Herzig, J.P., Leclerc, D.M. & Goff, P.L. 1970 Flow of suspensions through porous media–application to deep filtration. Ind. Engng Chem. Res. 62 (5), 8–35.CrossRefGoogle Scholar
Hestenes, M.R. 1969 Multiplier and gradient methods. J. Optim. Theor. Applics. 4 (5), 303–320.CrossRefGoogle Scholar
Hewitt, D.R., Daneshi, M., Balmforth, N.J. & Martinez, D.M. 2016 Obstructed and channelized viscoplastic flow in a Hele-Shaw cell. J. Fluid Mech. 790, 173–204.CrossRefGoogle Scholar
Hill, R. 1998 The Mathematical Theory of Plasticity, vol. 11. Oxford University Press.Google Scholar
Huilgol, R.R. 2015 Fluid Mechanics of Viscoplasticity. Springer.CrossRefGoogle Scholar
Jaisi, D.P., Saleh, N.B., Blake, R.E. & Elimelech, M. 2008 Transport of single-walled carbon nanotubes in porous media: filtration mechanisms and reversibility. Environ. Sci. Technol. 42 (22), 8317–8323.CrossRefGoogle ScholarPubMed
Kharabaf, H. & Yortsos, Y.C. 1997 Invasion percolation with memory. Phys. Rev. E 55 (6), 7177.CrossRefGoogle Scholar
Lester, D.R., Metcalfe, G. & Trefry, M.G. 2013 Is chaotic advection inherent to porous media flow? Phys. Rev. Lett. 111 (17), 174101.CrossRefGoogle ScholarPubMed
Liu, C., De Luca, A., Rosso, A. & Talon, L. 2019 Darcy's law for yield stress fluids. Phys. Rev. Lett. 122 (24), 245502.CrossRefGoogle ScholarPubMed
Mewis, J. & Wagner, N.J. 2012 Colloidal Suspension Rheology. Cambridge University Press.Google Scholar
Pandey, A. 2003 Solid-state fermentation. Biochem. Engng J. 13 (2–3), 81–84.CrossRefGoogle Scholar
Papanastasiou, T.C. 1987 Flows of materials with yield. J. Rheol. 31 (5), 385–404.CrossRefGoogle Scholar
Powell, M.J.D. 1978 Algorithms for nonlinear constraints that use Lagrangian functions. Math. Program. 14 (1), 224–248.CrossRefGoogle Scholar
Ronellenfitsch, H. & Dunkel, J. 2019 Chiral topological phases in designed mechanical networks. Front. Phys. 7, 178.CrossRefGoogle Scholar
Roquet, N. & Saramito, P. 2003 An adaptive finite element method for Bingham fluid flows around a cylinder. Comput. Meth. Appl. Mech. Engng 192 (31), 3317–3341.CrossRefGoogle Scholar
Roustaei, A., Chevalier, T., Talon, L. & Frigaard, I.A. 2016 Non-darcy effects in fracture flows of a yield stress fluid. J. Fluid Mech. 805, 222–261.CrossRefGoogle Scholar
Roustaei, A. & Frigaard, I.A. 2015 Residual drilling mud during conditioning of uneven boreholes in primary cementing. Part 2: steady laminar inertial flows. J. Non-Newtonian Fluid Mech. 226, 1–15.CrossRefGoogle Scholar
Sahimi, M. 2011 Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches. John Wiley & Sons.CrossRefGoogle Scholar
Saramito, P. 2007 A new constitutive equation for elastoviscoplastic fluid flows. J. Non-Newtonian Fluid Mech. 145 (1), 1–14.CrossRefGoogle Scholar
Saramito, P. 2016 Complex Fluids. Springer.CrossRefGoogle Scholar
Saramito, P. & Wachs, A. 2017 Progress in numerical simulation of yield stress fluid flows. Rheol. Acta 56 (3), 211–230.CrossRefGoogle Scholar
Schlesinger, W.H. 1999 Carbon sequestration in soils. Science 284 (5423), 2095.CrossRefGoogle Scholar
Silin, D. & Patzek, T. 2006 Pore space morphology analysis using maximal inscribed spheres. Phys. A 371 (2), 336–360.CrossRefGoogle Scholar
Sochi, T. 2005 Pore-scale modeling of non-Newtonian flow in porous media. PhD thesis, Imperial College London.Google Scholar
Stoop, N., Waisbord, N., Kantsler, V., Heinonen, V., Guasto, J.S. & Dunkel, J. 2019 Disorder-induced topological transition in porous media flow networks. J. Non-Newtonian Fluid Mech. 268, 66–74.CrossRefGoogle Scholar
Sun, H., Zhu, J., Baumann, D., Peng, L., Xu, Y., Shakir, I., Huang, Y. & Duan, X. 2019 Hierarchical 3D electrodes for electrochemical energy storage. Nat. Rev. Mater. 4 (1), 45–60.CrossRefGoogle Scholar
Talon, L., Auradou, H., Pessel, M. & Hansen, A. 2013 Geometry of optimal path hierarchies. Europhys. Lett. 103 (3), 30003.CrossRefGoogle Scholar
Talon, L. & Hansen, A. 2020 Effective rheology of bi-viscous non-Newtonian fluids in porous media. Front. Phys. 7, 225.CrossRefGoogle Scholar
Tien, C. & Payatakes, A.C. 1979 Advances in deep bed filtration. AIChE J. 25 (5), 737–759.CrossRefGoogle Scholar
Tjaden, B., Cooper, S.J., Brett, D.J.L., Kramer, D. & Shearing, P.R. 2016 On the origin and application of the bruggeman correlation for analysing transport phenomena in electrochemical systems. Curr. Opin. Chem. Engng 12, 44–51.CrossRefGoogle Scholar
Torquato, S. 1995 Nearest-neighbor statistics for packings of hard spheres and disks. Phys. Rev. E 51 (4), 3170.CrossRefGoogle ScholarPubMed
Torquato, S. & Haslach, H.W. Jr. 2002 Random heterogeneous materials: microstructure and macroscopic properties. Appl. Mech. Rev. 55 (4), B62–B63.CrossRefGoogle Scholar
Torquato, S., Lu, B. & Rubinstein, J. 1990 Nearest-neighbor distribution functions in many-body systems. Phys. Rev. A 41 (4), 2059.CrossRefGoogle ScholarPubMed
Torquato, S., Uche, O.U. & Stillinger, F.H. 2006 Random sequential addition of hard spheres in high Euclidean dimensions. Phys. Rev. E 74 (6), 061308.CrossRefGoogle ScholarPubMed
Tsamopoulos, J., Dimakopoulos, Y., Chatzidai, N., Karapetsas, G. & Pavlidis, M. 2008 Steady bubble rise and deformation in Newtonian and viscoplastic fluids and conditions for bubble entrapment. J. Fluid Mech. 601, 123–164.CrossRefGoogle Scholar
Waisbord, N., Stoop, N., Walkama, D.M., Dunkel, J. & Guasto, J.S. 2019 Anomalous percolation flow transition of yield stress fluids in porous media. Phys. Rev. Fluids 4 (6), 063303.CrossRefGoogle Scholar
Wei, T., Fan, F.Y., Helal, A., Smith, K.C., McKinley, G.H., Chiang, Y. & Lewis, J.A. 2015 Biphasic electrode suspensions for li-ion semi-solid flow cells with high energy density, fast charge transport, and low-dissipation flow. Adv. Energy Mater. 5 (15), 1500535.CrossRefGoogle Scholar
Zhang, G. & Torquato, S. 2013 Precise algorithm to generate random sequential addition of hard hyperspheres at saturation. Phys. Rev. E 88 (5), 053312.CrossRefGoogle ScholarPubMed
Zhu, Y., et al. . 2020 High-energy and high-power Zn-Ni flow batteries with semi-solid electrodes. Sustain. Energy Fuels 4, 4076.CrossRefGoogle Scholar
Figure 0

Figure 1.(a) Typical configuration of a porous medium of length $L$ and porosity $1-\phi$ that is made of non-overlapping monodisperse discs of radius $R_s$, and its network representation. Across the domain, a macroscopic pressure gradient ${\rm \Delta} P/L$ is applied. (b) Schematic of three network edges with different $h_i$ and $\ell _i (i=1,2,3)$. The network representation includes the local geometric characteristics of the complex porous medium structure. (c) Characteristic pore-size distribution for the network shown in (a) as derived from the network model. (d) Typical shear stress response $\tau$ as a function of the applied shear rate $\dot {\gamma }$ of a viscoplastic fluid with yield stress $\tau _y$.

Figure 1

Figure 2.Validation of the network model for the first open channel against simulations for a pressure-driven Bingham fluid in a porous medium with $R_s/L=0.02$ and $R_s/L=0.1$. (a) Simulation results of the open pathway for conditions near the critical pressure drop ${\rm \Delta} P_c$. The contour plot shows the magnitude of the local velocity, normalized with the maximum velocity across the channel. (b) Network model predictions for the first open channel. The cases of $\phi =0.3$ and $\phi =0.5$ are examined. It is clear that both the full Bingham fluid simulation and the network model predict the same location for the first open channel.

Figure 2

Table 1.Predicted normalized critical pressure drop ${\rm \Delta} P_c$ required to open the first open channel for the configurations shown in figure 2. We show both the predictions of the network model and the full simulation results, along with the relative error $({\rm \Delta} P_{c}^{{sim}}-{\rm \Delta} P_{c}^{{net}} ) / {\rm \Delta} P_{c}^{{sim}}$. The network model is adequate on predicting both the first open channel and the critical pressure drop required to open it.

Figure 3

Figure 3.Statistics of the predictions of the network model for 500 realizations per value of $\phi$. (a) Critical pressure drop ${\rm \Delta} P_c$ as a function of the volume fraction $\phi$ for different ratios of $R_s/L$; ${\rm \Delta} P_c$ is normalized with the yield stress of the fluid $\tau _y$. The error bars indicate the variance around the mean value. Inset – scaled critical pressure drop ${\rm \Delta} P_c R_s/\tau _y L$ in terms of $\phi$. For non-overlapping discs, the results for different $R_s/L$ collapse into a master curve. (b) Probability density distributions for the normalized arclength of the first open channel $L_c/L-1$. The cases of $R_s/L=0.02$ and $R_s/L=0.08$ are shown, respectively, for $\phi =0.1$ and $\phi =0.5$. In general, increasing $\phi$ leads to increase of ${\rm \Delta} P_c$ required to open the first channel, but also leads to a wide distribution of arclengths, which provide a large uncertainty in the critical macroscopic pressure gradient ${\rm \Delta} P_c/L_c$ required to yield the fluid in the porous medium.

Figure 4

Figure 4.Scalings of the geometric characteristics of a porous medium made of non-overlapping discs. (a) Inverse channel gap $\langle R_s/h \rangle$ thickness vs. $(1-\phi )^{-1}$. From our calculations it is clear that $\langle R_s/h \rangle$ scales with the inverse of the porosity $1-\phi$. (b) First open channel length $\langle L_c/L \rangle$ vs. $\phi$. Based on a simple estimate of the geometric tortuosity, we find that $\langle L_c/L \rangle$ scales linearly with $\phi$. (c) Non-dimensional pressure gradient ${\rm \Delta} \tilde {P}_c={\rm \Delta} P_c R_s/\tau _y L$ vs. $\phi /(1-\phi )$. According to (3.1), the product of $\langle L_c/L \rangle$ and $\langle R_s/h \rangle$ indicates that ${\rm \Delta} \tilde {P}_c$ scales with $\phi /(1-\phi )$.