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,
\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)
\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
\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
\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.
