Yielding to percolation: a universal scale

A theoretical and computational study analysing the initiation of yield-stress fluids percolation in porous media is presented. Yield-stress fluid flows through porous media are complicated due to the non-linear rheological behaviour of this type of fluids, rendering the conventional Darcy type approach invalid. A critical pressure gradient must be exceeded to commence the flow of a yield-stress fluid in a porous medium. As the first step in generalising the Darcy law for yield-stress fluids, a universal scale based on the variational formulation of the energy equation is derived for the critical pressure gradient which reduces to purely geometrical feature of the porous media. The presented scaling is then validated by both exhaustive numerical simulations (using an adaptive finite element approach based on the augmented Lagrangian method), and also the previously published data. The considered porous media is constructed by randomised obstacles with various topologies; namely, square, circular and alternatively polygonal obstacles which are mimicked based on Voronoi tessellation of circular cases. Moreover, computations for the bi-dispersed obstacle cases are performed which further demonstrate the validity of the proposed universal scaling.


Introduction
Yield-stress fluid flows through porous media are inherent to many industries including filtration, oil & gas and mining (Frigaard et al. 2017) and also numerous other applications such as biomedical treatments (Keating et al. 2003).Although in the case of Newtonian fluids many aspects of flows in porous media are well-discussed in the literature, when it comes to yield-stress fluids, our understanding of the phenomenon is limited mainly because modeling this problem is cumbersome due to the computational costs and/or the complexity of the experiments needed to carry out the analysis.
To overcome these barriers, several studies focused on pore-scale features of this problem (Bleyer & Coussot 2014;Shahsavari & McKinley 2016;De Vita et al. 2018;Bauer et al. 2019;Waisbord et al. 2019;Chaparian et al. 2020), however, it is yet unclear how to link/upscale the studies in micro-scale to macro-scale, especially due to the non-linearity of the constitutive equations which renders the bulk transport properties unpredictable from pore-scale dynamics.Nevertheless, in the intricate transport mechanism of yieldstress fluids through porous media, several mutual features can be identified regardless of the scales on which the previous studies are focused.In a number of studies (Talon & Bauer 2013;Liu et al. 2019;Chaparian & Tammisola 2021;Talon 2022), four regimes are detected in terms of flow rate (Q)-applied pressure gradient (∆P/L): (i) When the applied pressure gradient is less than a critical pressure gradient (∆P c /L), there is no flow (Q = 0); (ii) If the applied pressure gradient slightly exceeds the critical value, the flow is extremely localised in a channel and the flow rate linearly scales with the excessive pressure gradient where other parts of the fluid are quiescent; (iii) The third regime emerges when the applied pressure gradient increases, more and more channels appear (moderate values of pressure gradient) and the flow rate scales quadratically with the excessive applied pressure gradient; (iv) Finally when the applied pressure gradient is much higher than the critical value, the flow rate again scales linearly with the excessive pressure gradient.
Although these generic features/scales have been evidenced in a large number of studies, still the lack of an inclusive Darcy type expression for bulk properties is evident.The very first step for finding such a generic model is to thoroughly understand the pressure gradient threshold and more generally the yield limit which scaffolds any further progression of this aim.
In spite of the previous efforts to address the yield limit of the current problem (Liu et al. 2019;Chaparian et al. 2020;Fraggedakis et al. 2021), mostly the findings are case dependent, thereby limiting their application for more complicated practical systems.As discussed, in this limit the flow is extremely heterogeneous, hence, pore-scale studies are not fully reliable since they do not contain any statistical data in "real" porous media where a wider range of length scales are involved.Thus, in the present study, the aim is to derive a theoretical model based on yield-stress fluid flows principles and then validate the proposed model with exhaustive simulations.
To this end, we construct our porous media by randomly distributed obstacles of various shapes and lengths to avoid any biased results.Namely, three major types of obstacles are considered: circles, squares and polygons.Then, fluid flow simulations based on the adaptive augmented Lagrangian approach (Glowinski & Wachs 2011;Roquet & Saramito 2003) are performed which is shown to be a reliable tool for investigating the present problem, especially at the yield limit where non-regularised rheology is essential (Frigaard & Nouar 2005).To be fit for purpose, both mono-dispersed and bi-dispersed systems are considered.We have recently delved into mono-dispersed circular obstacles by the means of pore-network approaches where a large data set has been generated (Fraggedakis et al. 2021).This data set is adopted here in conjunction with the present computational data for further validation of the proposed theory.
The outline of the present paper is as follows: the problem is described in §2.1 and the details of the utilized numerical method and porous media construction are highlighted.The numerical results are depicted in §3.The theory is developed in §4 and the comparison with the computational results is performed.Conclusions are drawn in §5.

Mathematical formulation
We consider incompressible two-dimensional Stokes flow through a set of obstacles (i.e.X) in a box of size L × L (i.e.Ω) which is governed by, where p, τ, and u represent the pressure, deviatoric stress tensor, and the velocity vector of the fluid, respectively.We use the Bingham model to describe the fluid's rheology, in which γ is the rate of strain tensor (i.e.∇u + ∇u T ) and ∥ • ∥ is the second invariant of the tensor.Therefore, yielding obeys the von Mises criterion.
The above equations are non-dimensional and B = τy l/μ V is the Bingham number, where μ is the plastic viscosity of the Bingham fluid, V is the mean inlet velocity and l is the characteristic length scale which will be fixed later in §2.2.Hence, the Bingham number is the ratio of the yield stress of the fluid to the characteristic viscous stress.To derive the equations (2.1) and (2.2), we use the following scalings: where x and y are the coordinates in the streamwise and spanwise directions, respectively (see figure 1).Please note that all the variables with hat are dimensional throughout the paper; the same symbols are used for the dimensionless parameters without •.As mentioned above, V is the mean inlet velocity, hence, where Q is the flow rate and L inl is the length of the domain's inlet, i.e. the obstructed length by the solid obstacles is subtracted from L to calculate L inl (see figure 1).Therefore, in this setting, the flow rate is always equal to L inl , irrespective of the Bingham number.This approach in formulating the present problem is called Resistance formulation or [R].Indeed, the yield limit in this type of problem setup moves to B → ∞.
We will predominantly use this approach in our following simulations and analytical derivations.Alternatively, another formulation is possible: In [M] approach, the applied pressure gradient is used to scale the pressure and the stress tensor (i.e.∆ P L l), while the velocity vector is scaled with l2 μ ∆ P L .Hence, as a result, the non-dimensional applied pressure gradient in [M] is always equal to unity.
In [M] formulation, the independent flow parameter is, which is known as the yield number.Indeed, the flow rate changes as the yield number varies: it is zero when Y ⩾ Y c and increases as the yield number drops below Y c and decreases.Indeed, the yield limit in [M] is marked by Y c which is the critical yield number; if Y < Y c , the applied pressure gradient is enough to overcome the yield stress resistance and the fluid flows inside the medium.
There is a one-to-one map between [R] and [M] approaches: these two distinct formulations are linked together with Y (∆P/L) = B.This makes the interpretation of the results feasible; no matter the analysis (analytical, computational, etc.) is done in [R] or [M] settings.For more detailed explanations of these two formulations in porous media flows or more general pressure-driven flows, readers are refereed to Chaparian & Tammisola (2021).

Porous media construction
To construct the porous media for the fluid flow simulations, we randomly distribute non-overlapping obstacles (X) inside a square domain (Ω) of size L × L = 50 × 50; see figure 2. Indeed, the centre of each obstacle is chosen randomly with uniform distribution in the interval [−ϵ, L + ϵ] × [−ϵ, L + ϵ] and then it will be checked if the obstacle satisfies the non-overlapping condition.Here ϵ is introduced to let the obstacles cross the computational borders.
Three different obstacle topologies are used: circles, squares, and polygons.We consider mono-dispersed and bi-dispersed cases.In the mono-dispersed circular cases, the radius of the obstacles is used as the length scale, l = R, which deduces each individual obstacle area to be equal to π.In the mono-dispersed square cases, to be consistent, the individual area of an obstacle is again π or indeed l = Ls /π where Ls is the length of squares' edges.
In the bi-dispersed cases, the area of the larger obstacles is 25π while the area of the smaller ones is still π.This is the only parameter which is fixed in the construction of bi-dispersed cases.To ensure generality of the results, both the positions and the number of the larger obstacles are also chosen completely randomly.
For the case of polygons (see panels (c,f) in figure 2), firstly the domain is partitioned based on the Voronoi tessellation in which the centres of circles are adopted as the set of points in the Euclidean plane.Then each Voronoi cell (the edges of cells are depicted in red in figure 2) is squeezed (or expanded in the bi-dispersed cases) to get the desired area of π (or 25π in the bi-dispersed cases).Hence, this method provides us with a variety of shapes for the polygon cases.
As mentioned, here we are interested in 2D flows, hence, the solid "volume" fraction in the porous media is denoted by ϕ = meas(X)/meas(Ω).Therefore, the porosity of the medium (i.e.void fraction) can be represented simply by 1 − ϕ.
Note that in the polygon bi-dispersed cases, the obstacles may weakly overlap because of the expansion of the cells associated with the larger obstacles.In these cases, the effective solid volume fraction is considered.

Computational details
We implement augmented Lagrangian method to simulate the viscoplastic fluid flow (Glowinski & Wachs 2011;Roquet & Saramito 2003).This method is capable of handling the non-differentiable Bingham model by relaxing the rate of the strain tensor.An open source finite element environment-FreeFEM++ (Hecht 2012)-is used for discretisation and meshing which has been widely discussed and validated in our previous studies; for more details (choice of elements etc.) please see Chaparian & Frigaard (2017); Iglesias As discussed in §2.1, in [R] setting, the flow rate must be equal to L inl , hence, the imposed pressure gradient ∆P/L (which is a body force term in the numerical implementation) will be iterated to match the flow rate (Roustaei et al. 2015).
In the present study, a number of simulations are performed at different porosities to validate the scaling which will be derived in §4.As mentioned in §1, the main aim here is presenting a mathematical model for the yield limit based on physical features of the problem, and then validate it with the present simulations and the previous published data.Due to high computational cost of the full fluid flow simulations, we do not follow a statistical approach by simulating the flow in many realisations here.Rather, we use our data previously published in

Universal scale
For the present problem defined in §2.1, the energy equation at the steady state implies that the work done by the applied pressure gradient (i.e.(∆ P / L) Ω\ X û d Â) balances the total dissipation (i.e.Ω\ X (τ : γ) d Â = μ Ω\ X ( γ : γ) d Â + τy Ω\ X ∥ γ∥ d Â) which in dimensionless form reads, where a(u, u) is the viscous dissipation and B j(u) is the plastic dissipation.At the yield limit ), the viscous dissipation (which is quadratic in terms of γ) is at least one order of magnitude less than the plastic dissipation (Frigaard 2019;Chaparian et al. 2020), hence, the critical yield number (or indeed the inverse of the non-dimensional critical pressure gradient) can be predicted by, One can re-write the numerator as, since the flow rate is equal to L inl ; see expression (2.3).At the yield limit, the flow in the porous media is localised to a single channel.Thus, to find the scalings for j(u) and L inl at this limit, it is worth revisiting the two-dimensional Poiseuille flow of a yield-stress fluid.In this type of flow, the fluid moves as a core unyielded region with a constant velocity which is sandwiched between two sheared regions in which the velocity profile is parabolic.In the yield limit, these two sheared regions are viscoplastic boundary layers (Piau 2002;Balmforth et al. 2017) In our recent study (Fraggedakis et al. 2021), we have shown that the mean height of the first open channel scales with the porosity, i.e. ⟨h ch ⟩ ∼ 1 − ϕ and the mean relative length of the first channel scales with the volume fraction, i.e. ⟨L ch ⟩ /L ∼ ϕ, where ⟨•⟩ stands for the mean quantity which is acquired by ensemble averaging through different simulations and also various porosities.To elaborate, in a condensed system of obstacles (i.e.low porosities), h ch is smaller since the fluid path is squeezed between the obstacles or the mean void length between the obstacles becomes smaller as the solid volume fraction increases.On the other hand, the mean relative length of the first channel or tortuosity (i.e.L ch /L) scales with the solid volume fraction since in a denser system, the minimum path's shape is zigzag rather than straight which is more probable in a more dilute system of obstacles.These interpretations are evidenced in figure 5  Inserting the scales for the mean height and the mean relative length of the first channel to expression (4.3), the critical yield number can be re-written as: which means that the critical yield number scales with the ratio of the void space to the solid (i.e.obstructed) space.
In figure 6, we present a comparison of the theory (i.e.expression (4.4)) with the data associated with the simulations performed in the current study and also the previously published data: the non-dimensional critical pressure gradient (i.e.1/Y c ) is plotted versus ϕ/(1 − ϕ).The dashed orange line is the scale derived above, i.e. expression (4.4).The hollow symbols are the present computed data: black and purple colours are devoted to mono-dispersed and bi-dispersed cases, respectively.Circles, squares, and pentagrams represent the circle, square, and polygon obstacles, respectively.The filled circle symbols with the uncertainty bars are the data borrowed from Fraggedakis et al. (2021) where a pore-network approach is utilised to analyse a large number of realisations (∼ 500 for each porosity) with circular obstacles where each colour represents an specific R/ L ratio.Indeed, the filled circle symbols are the ensemble averages of all previously performed simulations and the uncertainty bars represent the range of obtained values.For more clarification of the used data, please see Fraggedakis et al. (2021).However, as explained in §2.3, the current data is acquired through individual simulations (i.e. they are not ensemble averages of many simulations), hence, no uncertainty bars are associated with the new data (i.e.hollow symbols).
A reasonable agreement can be observed between the derived scale (with a fitted slope ≈ 3.14 or π) and the computational data for all class of considered topologies.Moreover, the bi-dispersed cases data also fits reasonably to the proposed theory.
In a very recent study, using "variational linear comparison" homogenisation method, Castañeda (2023) has derived an upper-bound for the critical pressure gradient where the solution of Newtonian fluids used as a test function in the dissipation-rate potential of viscoplastic fluids.This upper-bound is shown by the cyan line in the inset of figure 6 along with the proposed universal scale for comparison.Note that the upper-bound proposed by Castañeda (2023) has a linear functionality with ϕ/(1 − ϕ) which further validates the universal scale derived here, although its slope is steeper which is not surprising as it is an upper-bound.(2021).Each colour intensity is dedicated to a different value of R/ L between 0.02 to 0.1 (see the reference for more details).The black and purple hollow symbols devote the mono-dispersed and bi-dispersed cases, respectively.Circles, squares, and pentagrams represent the circle, square, and polygon obstacles, respectively.Inset: comparison between the upper-bound of the critical pressure gradient (cyan line) derived by Castañeda (2023) and the proposed universal scale (dashed orange line).Please note that the axes of the inset are the same as the main figure.

Concluding remarks
Adaptive finite element simulations based on augmented Lagrangian scheme were performed to study the fluid flows of yield-stress fluids in porous media.The specific objective was to fully understand the yield limit of this type of flows and propose a theory to address the critical applied pressure gradient which should be exceeded for flow assurance purposes.This is a vital and a very base step in proposing a generic Darcy type expression for bulk transport properties of the yield-stress fluid flows in porous media.
For this aim, and to avoid prevailing analysis, flows in various porous media constructed with a wide range of obstacle shapes are investigated.The studied geometries have been generated by randomly distributing non-overlapping obstacles of circular and square shapes.In addition, more complicated topologies (i.e.polygon obstacles), have been generated by using the Voronoi tessellation of circular cases.The computational data includes both mono-dispersed and bi-dispersed systems.
In the yield limit, which is the main focus of the present study, the flow is restricted to a single channel connecting the inlet to outlet, while the fluid outside of it is unyielded and thus quiescent.The configuration of this very first channel has been investigated in our previous study (Fraggedakis et al. 2021) and statistical geometrical properties (e.g.height and length) are reported as a function of the solid volume fraction (ϕ) or alternatively the porosity of the domain (1−ϕ) which can be summarised as ⟨h ch ⟩ ∼ 1−ϕ and ⟨L ch ⟩ /L ∼ ϕ.
A theory was proposed based on variational formulation of the energy equation.
The leading order plastic dissipation has been approximated by a channel Poiseuille flow at the yield limit where the channel dimensions are borrowed from the discussed statistical results (Fraggedakis et al. 2021).Indeed, in the very first channel, the transport mechanism is predominantly postulated by the core unyielded plug in the middle of the channel and the leading order plastic dissipation occuring in the sheared boundary layer between the quiescent fluid outside of the channel and the mobilised core unyielded region.It should be noted that due to the complex shape of this limiting channel in the porous media, the mechanism is not as simple as explained above since the limiting channel is not straight and channel height varies (especially in the dense systems); see figure 5. Thus, the core unyielded plug and the adjacent boundary layers are not uniform.Nevertheless, since the mean height and length of the channel is used in our model, the proposed scaling is still valid in the leading order.This has been assessed using the obtained computational data for a wide range of obstacle topologies mentioned above and also previously published data.We have shown that our theoretical approach is capable of predicting the numerical data with a reasonable agreement.Due to the high cost of unregularised numerical simulations of yield-stress fluid flows and also handling various shapes of obstacles, the available data, especially in the yield limit, is limited.This limitation is more evident in three-dimensional flows.Although in some studies (Bittleston et al. 2002;Pelipenko & Frigaard 2004;Hewitt et al. 2016;Izadi et al. 2023), the Hele-Shaw approximation for yield-stress fluids has been developed, still the lack of a compelling study linking this pore-scale approximation to bulk transport mechanisms/features in 3D is evident.This is left for future investigations, both theoretically and computationally, which is a massive step forward for many industrial applications.

Figure 1 .
Figure 1.Schematic of the coordinate system directions and the inlet length L inl which in this case consists of two segments depicted in blue.

Figure 2 .
Figure 2. Schematic of the porous media: top row are mono-dispersed topologies and bottom row are bi-dispersed ones.(a,d) square obstacles; (b,e) circular obstacles; (c,f) generated polygon obstacles based on Voronoi tessellation of panels (b,e).
Fraggedakis et al. (2021) which will be discussed later in §4.Recent advances in computational methods of viscoplastic fluids (e.g.known as PAL & FISTA methods) accelerate the simulations of this type of fluids, yet the implementation of these methods is beyond the scope of this work.Interested readers are refereed to Dimakopoulos et al. (2018) and Treskatis et al. (2016, 2018).

Figure 3 .
Figure 3. Mesh generation for a sample case: (a) initial mesh ("uniform" coarse grid), (b) final mesh after 6 cycles of adaptation.This mesh is associated with the simulation illustrated in panel (d) of figure 4. Note that only part of the mesh in the white window of panel (d) of figure 4 (at the pore-scale) is shown here.

Figure 4
Figure 4 shows the flow in the six sample geometries at ϕ = 0.45 & B = 10 3 .As discussed in §2.1, the yield limit in the [R] setting goes to B → ∞, so in the illustrated examples for this relatively large Bingham number, the channelisation is clear.However, clearly for different geometries, different "large" Bingham numbers are required to get only the very first open channel.This translates to different critical yield numbers which is expected for different topologies and will be discussed in §4 with other features of the flows.

Figure 4 .
Figure 4. Contour of velocity (i.e.|u|) for 6 sample simulations at ϕ = 0.45 & B = 10 3 .Top row panels are mono-dispersed cases and the bottom panels are the bi-dispersed ones.The white window in panel (d) marks where the mesh represented in figure 3 belongs to.
with thickness δ.To simplify the plastic dissipation functional j(u) substantially, we approximate the flow in the first open channel with the discussed Poiseuille flow.Hence, the leading order of ∥ γ∥ can be approximated as ≈ 2(U ch /δ) δ L ch ∼ U ch L ch in the boundary layers where the index ch stands for the first open channel.Indeed, U ch and L ch represent the core unyielded region velocity and the length of the first channel, respectively; see figure5(a).Moreover, the continuity equation in the leading order obeys Q = L inl ≈ U ch h ch which allows us to rewrite U ch in terms of the flow rate and the channel height.Thus,

Figure 6 .
Figure 6.Comparison between our theory and the computational result: non-dimensional critical pressure gradient versus ϕ/(1 − ϕ).The dashed orange line is the scale derived in (4.1).The filled circle symbols with uncertainty bars are the data borrowed from Fraggedakis et al.(2021).Each colour intensity is dedicated to a different value of R/ L between 0.02 to 0.1 (see the reference for more details).The black and purple hollow symbols devote the mono-dispersed and bi-dispersed cases, respectively.Circles, squares, and pentagrams represent the circle, square, and polygon obstacles, respectively.Inset: comparison between the upper-bound of the critical pressure gradient (cyan line) derived byCastañeda (2023) and the proposed universal scale (dashed orange line).Please note that the axes of the inset are the same as the main figure.