Numerical study of Bingham flow in macrosopic two dimensional heterogenous porous media

TL;DR

This study numerically investigates the flow behavior of yield stress non-Newtonian fluids in two-dimensional heterogeneous porous media, revealing fractal flow paths and scale-dependent properties.

Contribution

It provides new insights into the macroscopic flow patterns of yield stress fluids in heterogeneous media, highlighting scale-dependent fractal characteristics and independence from heterogeneity amplitude.

Findings

Flow paths exhibit fractal features with scale exponents.

Exponents are independent of heterogeneity amplitude for log-normal distributions.

Macroscopic exponents differ from microscopic ones, indicating different governing equations.

Abstract

The flow of non-Newtonian fluids is ubiquitous in many applications in the geological and industrial context. We focus here on yield stress fluids (YSF), i.e. a material that requires minimal stress to flow. We study numerically the flow of yield stress fluids in 2D porous media on a macroscopic scale in the presence of local heterogeneities. As with the microscopic problem, heterogeneities are of crucial importance because some regions will flow more easily than others. As a result, the flow is characterized by preferential flow paths with fractal features. These fractal properties are characterized by different scale exponents that will be determined and analyzed. One of the salient features of these results is that these exponents seem to be independent of the amplitude of heterogeneities for a log-normal distribution. In addition, these exponents appear to differ from those at the microscopic level, illustrating the fact that, although similar, the two scales are governed by different sets of equations.

Equations

$$\dot{\gamma}=\left\{\begin{array}[]{lcc}0&\rm{if}&\tau<\tau_{0}\\ (\frac{1}{C}(\tau-\tau_{0}))^{1/n}&\rm{if}&\tau>\tau_{0}\\ \end{array}\right.,$$

$$\vec{u}=\left\{\begin{array}[]{lcc}-\frac{k}{\eta_{0}}({\color[rgb]{0,0,0}\vec{\nabla}}P-{\color[rgb]{0,0,0}G}\frac{{\color[rgb]{0,0,0}\vec{\nabla}}P}{|{\color[rgb]{0,0,0}\vec{\nabla}}P|})=-\frac{k}{\eta_{0}}({\color[rgb]{0,0,0}\vec{\nabla}}P+{\color[rgb]{0,0,0}G}\frac{\vec{u}}{|\vec{u}|})&if&|{\color[rgb]{0,0,0}\vec{\nabla}}P|>{\color[rgb]{0,0,0}G}\\ \vec{0}&if&|{\color[rgb]{0,0,0}\vec{\nabla}}P|<{\color[rgb]{0,0,0}G}\end{array}\right.,$$

$$\vec{0}=-{\color[rgb]{0,0,0}\vec{\nabla}}P-\frac{\eta_{0}}{k}\left(1+\frac{k{\color[rgb]{0,0,0}G}}{\eta_{0}|u|}\right)\vec{u}.$$

$$\vec{0}=-{\color[rgb]{0,0,0}\vec{\nabla}}P-\frac{{\color[rgb]{0,0,0}\eta_{eff}}(|u|)}{k}\vec{u},$$

$$\vec{0}=-{\color[rgb]{0,0,0}\vec{\nabla}}P-\frac{{\color[rgb]{0,0,0}\eta_{eff}}(|u|)}{k}\vec{u}+\eta_{B}\Delta\vec{u},$$

$$pdf(f=\log k)=\frac{1}{\sqrt{2\pi \sigma }}\exp -\frac{(f-f_0{)}^2}{2{\sigma }^2},$$

$$\int f(\vec{r}+{\vec{r}}_0)f({\vec{r}}_0))dxdy=e^{-\frac{\pi }{8{\lambda }^2}(x^2+y^2)}$$

$$p(\ln{{\color[rgb]{0,0,0}G}})=\sqrt{\frac{2}{\pi\sigma}}\exp{(-\frac{2(\ln{{\color[rgb]{0,0,0}G}}-\ln{A/\sqrt{k_{0}}})^{2}}{\sigma^{2}})}.$$

$$X=\lambda \tilde{X},P=\lambda A/\sqrt{k_0}\tilde{P}\mathrm{and}\mathrm{u}=\frac{\mathrm{A}\sqrt{{\mathrm{k}}_0}}{{\eta }_0}\tilde{\mathrm{u}},$$

$$f_q(\vec{r}+{\vec{c}}_q,t+1)={\hat{f}}_q(\vec{r},t),q=0\dots 8,$$

$${\hat{f}}_q(\vec{r},t)=[f_q-s^{+}{n}_q^{+}-s^{-}{n}_q^{-}],q=0\dots 4$$

$${\hat{f}}_{\bar{q}}(\vec{r},t)=[f_{\bar{q}}-s^{+}{n}_q^{+}+s^{-}{n}_q^{-}],q=0\dots 4,$$

$$n_q^{\pm }$$

$$P=\rho c_s^2\mathrm{with}\rho =\sum _0^8{\mathrm{f}}_{\mathrm{q}}$$

$$\vec{u}=\frac{2\vec{J}}{2+\frac{{\color[rgb]{0,0,0}\eta_{eff}}(\vec{u})}{k}}\;\;\rm{with}\;\;\vec{J}=\sum_{0}^{8}f_{q}\vec{c}_{q},$$

$$e_q^{+}=t_{q}P,e_0=\rho -\sum _1^8{e}_q^{+}$$

$$e_q^{-}=t_q(\vec{u}.{\vec{c}}_q+{\Lambda }^{-}{F}_q).$$

$$F_{q}=-t_{q}\frac{{\color[rgb]{0,0,0}\eta_{eff}}(\vec{u})}{k}\vec{j}.\vec{c}_{q}.$$

$$s^{+}=\frac{2}{2\tilde{\Lambda}^{+}+1}\;\;\rm{with}\;\;\tilde{\Lambda}^{+}=\frac{9(4+\frac{{\color[rgb]{0,0,0}\eta_{eff}}(\vec{u})}{\eta_{0}k})}{4(3+2\Lambda\frac{{\color[rgb]{0,0,0}\eta_{eff}}(\vec{u})}{\eta_{0}k})}\eta_{0}$$

$$s^{-}=\frac{2}{2{\Lambda }^{-}+1}.$$

$$\kappa (x,y)=\sin (2\pi y)\frac{{\kappa }_{max}-{\kappa }_{min}}{2}+2{\kappa }_{min}.$$

$$Err=\frac{Q_{LBM}-Q_{FD}}{Q_{FD}}.$$

$${\color[rgb]{0,0,0}\Delta\tilde{P_{c}}}=\min_{\cal C}\int{\color[rgb]{0,0,0}\tilde{{\color[rgb]{0,0,0}G}}}(\vec{r})ds,$$

$$\tilde{q}=B(\Delta\tilde{P}-{\color[rgb]{0,0,0}\Delta\tilde{P_{c}}})^{\color[rgb]{0,0,0}\alpha}.$$

$$P(\tilde{S})\propto\tilde{S}^{-{\color[rgb]{0,0,0}\tau}}e^{-\tilde{S}/{\color[rgb]{0,0,0}\tilde{S_{0}}}},$$

$${\color[rgb]{0,0,0}\tilde{S_{0}}}\propto\tilde{q}^{-{\color[rgb]{0,0,0}\gamma}}.$$

$$Bx^{\color[rgb]{0,0,0}\zeta}(1-x)^{{\color[rgb]{0,0,0}\zeta}},$$

$${\color[rgb]{0,0,0}\tilde{L}}\propto{\color[rgb]{0,0,0}\tilde{W}}^{\color[rgb]{0,0,0}\zeta},$$

Full text

1 Introduction

The flow of non-Newtonian fluids is ubiquitous in many applications in geological and industrial context. In this study, we focus on yield stress fluid (YSF), viz. material that requires a minimum amount of stress to flow. When the stress is below a critical value, the material behave as solid. Once the stress is above this threshold, the material can be sheared and thus flow.

These materials have many implications for industrial processes. Indeed mud, polymers, oil, foam suspension may present a yielding threshold [1]. The most important is certainly the extraction of heavy oil and has been the subject of many studies since the 1960s [2]. Yield stress fluids are also used for enhanced oil recovery where foam or polymers are injected into the ground to block preferential flow paths [3, 4]. Another application is the hydraulic fracking where the yield stress properties is used to prevent the closing of the fractures [5].

Because of the many applications, flow of yield stress fluids has been investigated in many studies [2, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Most of these studies focused on flow at the microscopic (pore) scale. The main objective was to establish a constitutive equation to relate the flow within the medium to the applied pressure drop, which represents a generalization of the Darcy’s law to yield stress fluids.

Yield stress fluids are usually described by the Herschel-Bulkley model which relates the shear rate, $\dot{\gamma}$, to the applied shear stress $\tau$:

[TABLE]

where $\tau_{0}$ is the yield stress, $C$ the consistency and $n$ an exponent. Here, we choose to restrict to the Bingham rheology which corresponds to $n=1$ (Bingham fluid) and $C=\eta_{0}$, the dynamic viscosity.

The usual governing equation for the flow of yield stress fluid at the field scale has been proposed by Entov [2] and Pascal [15]. They proposed a modified Darcy’s equation (see also [16, 17]) in the form:

[TABLE]

where $k$ is the permeability and {\color[rgb]{0,0,0}G} a positive scalar representing the limiting pressure gradient related to the yield stress. Here, we have used the fact that the velocity vector and the gradient of pressure are opposed direction. Since its introduction, this law has been validated by several experimental studies (see for instance [6, 7, 10]). It is important to note that both the limiting pressure and the permeability depend on the local topology of the material. These quantities may thus be subject to spacial variations at the macroscopic geological scale.

In recent studies [11, 18, 19], we have, however, demonstrated that, if the medium in sufficiently disordered, the flow curve might be non-linear close to the limiting pressure. Because of the structural disorder, some channel paths may be easier to flow than others. For a certain range of pressure, there is thus a non-linear increase of flowing paths responsible for the non-linearity of the flow rate curve. At sufficiently high pressure, once all the fluid has yielded, the flow curve becomes linear again and eq. (2) is valid.

One of the characteristics of these study was to demonstrate that the flow curve and structure were governed by a power law. Indeed, the flow rate satisfies Q\propto(\nabla P-{\color[rgb]{0,0,0}G})^{\alpha}, where $\alpha$ is an exponent close to $2$. Moreover, the size ($S$) distribution of the cluster of the fluid at rest follow : $P(S)\propto S^{-\tau}\exp{S/S_{0}}$ with $S_{0}\propto Q^{-\gamma}$. This distribution law characterizes the multi-scale aspect (fractal) of this problem. It is also reminiscent to other problems like percolation [20] or the avalanche dynamic [21, 19, 22]. It should also be noted that the exponent $\alpha$, $\tau$ and $\gamma$ were found to be independent of the type of disorder.

The aim of this study is to investigate the flow of a Bingham fluid at a macroscopic level, i.e by solving a modified Darcy’s law in a heterogenous permeability field. As noted above, the permeability and the limiting pressure gradient depends on the local topology of the medium. We therefore investigate the effect of spatial variations in these quantities on the flow. For example, one could expect some similarities with the microscale aspect such as the appearance of preferential paths. In their interesting work [23], Hewitt et al. studied Bingham fluids in the Hele-Shaw cell in the presence of obstacles where they observed an increase in preferential paths with pressure in fractures (Hele-Shaw with open variations).

The article is divided as follows. First, we will present the problem, the flow equations and the description of the porous medium. Secondly, we will present the numerical methods based on a Lattice Boltzmann scheme with two relaxation times. We will then present and discuss the results and the conclusion.

2 Governing equations

2.1 Darcy’s law

The purpose of this paper is to study the flow of Bingham fluid at a macroscopic scale. To this goal, we solve the modified Darcy’s equation which relates the local mean flow rate to the local mean gradient of pressure.

In the introduction, we have presented the standard Darcy’s equation (2) proposed in the literature. As discussed, this equation is not true close to the limiting pressure gradient, where one should expect a transitory regime which is more complex. However, because the transition is not yet fully understood and because it appears in a short range of pressure, for the sake of simplicity, we will assume that the standart eq. (2) is valid.

This equation can then be reformulated homogeneous to a balance of momentum:

[TABLE]

Here, we can see that the contribution of the yield stress can be seen as a constant body force opposed to the flow.

We have thus,

[TABLE]

where {\color[rgb]{0,0,0}\eta_{eff}}(|u|)=\eta_{0}\left(1+\frac{k{\color[rgb]{0,0,0}G}}{\eta_{0}|u|}\right) is a scalar effective viscosity that depends on the local velocity field. The model of effective viscosity is also very common to solve non-linear rheology in porous media. Following the experimental work of Hirasaki and Pope [24] or the one of Chauveteau [25], the idea is to introduce an effective shear rate related to the local mean velocity. The usual Darcy’s law is then modified with an effective viscosity varying with the local mean velocity.

However, we have introduced two modifications to this equation. First, like the Stokes equation for Bingham fluids, the effective viscosity has a zero velocity divergence point which can lead to numerical instabilities. We have tackled this problem by capping the viscosity by a maximum value $\eta_{eff}^{max}$ which is several order of magnitude larger than $\eta_{0}$ (typically $10^{9}$ larger). Another common technique to overcome this problem would have been to regularize this function around zero with an exponential term as proposed by Papanastasiou [26]. We tried both methods which gave similar results but the second was much slower.

The second comment is a general problem of the Darcy equation, which is a first order differential equation, in disorder media. Indeed, any discontinuity in the permeability field then leads to a discontinuity in the velocity field, which is not physical but also introduces numerical instabilites. To overcome this problem, one possibility is to introduce a diffusive term into the momentum equation as proposed by Brinkman [27]:

[TABLE]

where $\eta_{B}$ is a coefficient that depends on the viscosity and the geometry. Here, for the sake of simplicity, we keep it constant with $\eta_{B}=\eta_{0}$.

2.2 Porous medium

For the heterogenous permeability field, we use a log-normal distribution that has been widely used since the work of Gehlar [28] or Dagan [29] to describe heterogeneity at macroscopic scale. The permeabilty map is distributed according to:

[TABLE]

where $f_{0}$ and $\sigma$ are the mean and the standard deviation of $\log(K)$.

The field $f=\log(k)$ is correlated with a Gaussian function, set by the correlation length $\lambda$:

[TABLE]

The field is generated using a standard Fast Fourier Method (see Yiotis et al. [30] for details). The permeability field $k(\vec{r})$ is thus parametrized by three parameters, the mean $f_{0}$, the amplitude of the heterogeneities $\sigma$ and the correlation length $\lambda$. Moreover, we define $k_{0}=e^{f_{0}}$

From a phenomenological aspect, local critical pressure is expected to be related to permeability. However, no clear relationship has been proposed in the literature. We have therefore chosen to relate them using dimensional arguments. Indeed, one expects that the critical pressure can be written as {\color[rgb]{0,0,0}G}=\tau_{0}/d, where $d$ is a geometrical factor that has the dimension of a length (the typical pore size) [31]. Since permeability has the dimension of a length to the power $2$, a natural relationship is to assume {\color[rgb]{0,0,0}G}=A/\sqrt{k}, where $A$ is a prefactor. This relationship is, of course, purely phenomenological and cannot be true for any types of porous media. We expect to hold for media with same topology such has mono disperse bead packing.

From the distribution of permeability eq. (6), the probability density function of the limiting pressure also follows a log-normal distribution:

[TABLE]

2.3 Non-dimensionalization

In the following, we will use dimensionless quantities. For this, we nondimensionalize any lengths with the correlation length, $\lambda$, any pressure with the average critical pressure gradient $A/\sqrt{k_{0}}\lambda$ and any velocities with the average permeability and the average critical pressure gradient :

[TABLE]

where $\tilde{.}$ indicates dimensionless quantities.

3 Numerical method

To solve the Brinkman equation eq. (5), we used an improved two-relaxation time Lattice Boltzman method (IBF-TRT) proposed by I. Ginzburg [32, 33]. The basic idea of the Lattice Boltzmann method is to discretize the particle velocity distribution function on a grid. Here we used a 2-dimensional scheme with 9 different velocities ($\vec{c}_{q}$). We then introduce the population ${f_{q}}$ as the density of particles moving with the velocity ${\vec{c}_{q}}$. The algorithm consists of a succession of two main steps. The first is the propagation step Eq. (9), where we move the density on the grid according to its velocity. The second is the collision step Eq. (10), where populations meeting at the same node are redistributed using a collision operator. This collision part is described as a relaxation toward an equilibrium state that depends on local macroscopic quantities (pressure, speed…).

The main idea of the TRT scheme is to decompose each population into a symmetrical and an antisymmetrical part. Each component relax toward their equilibrium with its own relaxation rate. We define the even and odd part respectively as $f_{q}^{+}=\frac{f_{q}+f_{\bar{q}}}{2}$ and $f_{q}^{-}=\frac{f_{q}-f_{\bar{q}}}{2}$. Denoting by $\bar{q}$ the direction opposite to $q$ ($\vec{c}_{q}=-\vec{c}_{\bar{\textbullet}q}$) and assuming $\vec{c}_{0}=\vec{0}$ and that the $8$ nonzero velocities are ordered so that the first four directions are opposite to the last four.

The propagation step is described as:

[TABLE]

and the collision step:

[TABLE]

with

[TABLE]

where $e_{q}^{\pm}$ are the equilibrium distributions. The model has several parameters, $c_{s}$ the numerical sound speed, $\nu_{0}$ the Brinkman viscosity and $\Lambda$ a numerical parameter, from which we define $\Lambda^{+}=3\eta_{0}$ and $\Lambda^{-}=\Lambda/(3\eta_{0})$. The determination of the equilibrium function requires to compute the local pressure $P$ and momentum $\vec{u}$:

[TABLE]

and

[TABLE]

where the effective viscosity has be determined using the local velocity at the previous step.

The equilibrium function are then defined as:

[TABLE]

and

[TABLE]

$F_{q}$ is the body force which takes into accound the Darcy drag force (second term in eq. (5)) and is defined by:

[TABLE]

Finally, the two relaxation rates are defined as:

[TABLE]

and

[TABLE]

The main difference between the standard BF-TRT and the improved IBF-TRT scheme is the use of $\tilde{\Lambda^{+}}$ instead of $\Lambda^{+}$ for the relaxation parameter $s^{+}$. Indeed, the standard BF-TRT scheme leads to an error in the viscosity which depends on the local parameter. The change of $\Lambda^{+}$ is to compensate for this error [32, 33].

In this study, the typical porous medium size is $1024\times 4096$ nodes. The typical duration time is between $2$ and $100$ hours with $32$ CPUs.

3.1 Validation

To validate our numerical code, we simulated a Bingham fluid in a stratified porous medium with permeability distributed according to a sinusoidal function in the normal direction ($y$) and invariant in the flow direction ($x$), namely:

[TABLE]

We used different numbers of grid nodes, $N_{y}\in\{5,8,16,32,64,128\}$, to vary the spatial resolution. To perform a quantitative analysis, we solve eq. (5) using a standard second order finite difference method discretized with a fixed finer resolution of $512$ points.

In Fig. 1, we have plotted the velocity profile for different applied pressure field for a fixed resolution $N_{y}=16$. As we can see in this figure, the difference with the finite difference is quite satisfactory. We also compute in Fig. 2 the error difference between the total flow rate $Q_{LBM}$ and $Q_{FD}$ from the Lattice Boltzmann and the finite difference method respectively:

[TABLE]

The error remains quite low. The maximum of errors is reached close to the critical pressure where both solutions tend to zero.

Considering eq. (7) with $\lambda=5$, we can estimate that our simulations correspond to a periodicity of $N_{y}\simeq 20$. We expect thus an error of $1-2\%$ close to the first path opening, and of less than $0.1\%$ above.

4 Results

4.1 Problem description

We simulate the flow in domain size of $4096\times 1024$. The permeability distribution is log-normal with a mean $k_{0}=e^{-1}$, a correlation length \color[rgb]{0,0,0}\lambda=5 and a mean square deviation {\color[rgb]{0,0,0}\sigma}^{2} ranging from $0.1$ to $3$. The constant $A$ was set to $A=10^{-5}$. A constant pressure difference $\Delta P$ is then imposed at the inlet and outlet of the domain. After the simulation has converged, we compute the mean flow rate: $\tilde{q}=\frac{1}{V}\int udV$, where $V$ is the domain size.

When the imposed pressure is too low, the whole domain is in the solid state which results in a zero mean velocity (very small in fact due to the capping of the viscosity). Above this macroscopic threshold, the fluid then starts to flow in few channels. Figure 3 displays snapshots of the flow field for different {\color[rgb]{0,0,0}\sigma} and different applied pressure above this threshold. As it can be seen, close to this threshold the number of flowing paths is small. In principle, it should reduce to a single path but, due to numerical precision, the precise critical point is difficult to determine, particularly for low {\color[rgb]{0,0,0}\sigma}. As the pressure is increased, the number of flowing paths increase with the applied pressure. We also observe in these snapshots that the flowing paths seem to be more tortuous as the heterogeneity parameter, $\sigma$, is increased.

In the following, we will first present the evolution of the mean flow characteristics by increasing the pressure. Then we will study some geometrical aspect of the flow field.

4.2 Onset of flow

The first characteristic is the critical global pressure difference required to initiate the flow. Because of the yield stress property, we observe a minimal pressure drop {\color[rgb]{0,0,0}\Delta\tilde{P_{c}}} below which there is no flow (within the regularization error bar). As already described in other papers [34, 11, 23], this pressure corresponds to the appearance of the first flowing path. It can be formally defined as being the path which has the minimal critical pressure:

[TABLE]

where the path $\cal C$ are taken among all possible paths connecting the inlet to the outlet.

Several comments must be made on this critical path. We expect that the path and the associated critical pressure should depend on the disorder since they result from the competition of two opposite effects.

On the one hand, the minimal path tends to connect regions of low {\color[rgb]{0,0,0}\tilde{{\color[rgb]{0,0,0}G}}}(\vec{r}) (i.e. high permeability). But on the other hand, as indicated by eq. (18), the total length of the path has a contribution. The optimal path results thus from a balance between the search of the most permeable areas and the cost of increasing the length of the path.

This competion can clearly be observed in Fig. 3, where we have plot the first path for different $\sigma$. At lower $\sigma$, the path is almost straight because the lower {\color[rgb]{0,0,0}\tilde{{\color[rgb]{0,0,0}G}}}(\vec{r}) regions are not low enough to compensate for an increase of length. As $\sigma$ increases, the optimal path becomes more tortuous because the value of the lowest {\color[rgb]{0,0,0}\tilde{{\color[rgb]{0,0,0}G}}}(\vec{r}) regions decreases, which could be worth the detour.

We may note that the tortuosity of the first path may not necessarily increase with the amplitude of the disorder. A simple example would be to multiply the {\color[rgb]{0,0,0}\tilde{{\color[rgb]{0,0,0}G}}} field by a constant value (or by changing $A$), which has the effect of increasing the standard deviation but not the shape of the optimal path (see eq. (18)).

Finally, we should also note that the path selection of eq. (18) is very close to a standart problem in statistical physics named directed polymer problem [35], which consist in the finding of a minimal path in an energy lanscape. The minimization is however perform among ”directed” paths, meaning that their slope is bounded to avoid any overhangs. In this context, the directed polymer is known to be self-affine with a roughness exponent equal to $2/3$. This exponent is ”universal” in the sense that it does not depend on the distribution of the energy (see [35, 36]), provided that the distribution is not too ”extreme” (like power law, fat tail, distributions [37]). The roughness of our flowing paths will be discussed in section 4.4.

In Fig. 4, we plot the value of {\color[rgb]{0,0,0}\Delta\tilde{P_{c}}} as function of the amplitude of the disorder. As expected, {\color[rgb]{0,0,0}\Delta\tilde{P_{c}}} decreases with the amplitude of the disorder, starting from {\color[rgb]{0,0,0}\Delta\tilde{P_{c}}}(\sigma=0)={\color[rgb]{0,0,0}\tilde{{\color[rgb]{0,0,0}G}}}{\color[rgb]{0,0,0}\tilde{L}}, where {\color[rgb]{0,0,0}\tilde{L}} is the dimensionless total length of the system. Moreover, as indicated in the figure, the critial pressure seems to follow a gaussian function with {\color[rgb]{0,0,0}\sigma}.

Note that for low value of $\sigma$, the critical pressure and the path associated to it become more difficult to determine precisely. Indeed, since the disorder is small, all the possible paths are quite similar in term of pressure threshold. As a result, all flow paths appear very quickly in a narrow pressure range.

4.3 Flow regimes

We now study the evolution of the flow rate by increasing the pressure difference above the threshold. We have plotted in Fig. 5, the mean flow rate $q$ as function of the applied pressure substracted by the critical pressure \Delta\tilde{P}-{\color[rgb]{0,0,0}\Delta\tilde{P_{c}}}. As it can be seen, the mean debit follows a power-law over a certain range of pressure:

[TABLE]

More remarkable, we note that the exponent seems to be constant with the heterogenity paramter {\color[rgb]{0,0,0}\sigma} (within the error bar) {\color[rgb]{0,0,0}\alpha}=2.8\pm 0.05.

This flow behaviour is reminiscent of what has been observed at the pore scale [18]. However, there is a major difference that lies in the value of the exponent which is significantly higher ({\color[rgb]{0,0,0}\alpha}\simeq 2 at the pore scale). We believe that this is due to the contribution of the increase in the number of flow paths but also to the increase of the flow rate in each of them. As in the case of pore scale, a significant contribution is due to the increase in the number of flow paths with the applied pressure. However, while at the pore scale, the flow rate of each path is expected to increase linearly with pressure, this is not necessarily the case at the macroscopic level because the paths can also widen. Indeed, since the permeability field and the limit pressure are correlated, it is also expected that regions near open paths will also be easier to flow. The width of each individual path therefore increases with pressure, similarly to the stratified permeability distribution case considered in the validation section. This effect was quantified in Fig. 6 where we measured the distribution of the channel width $\tilde{W}_{c}$ as a function of the applied pressure. We can clearly see the broadening of the distribution with the increase in pressure.

Finally, for large enough pressure, once most of the domain is flowing, the number of the paths and their width cannot increase anymore, the flow rate then recover a linear behavior \tilde{q}\propto(\Delta\tilde{P}-{\color[rgb]{0,0,0}\Delta\tilde{P_{c}}}).

It should be noted that the amplitude of the disorder modifies the pressure range over which the intermediate regime is observed. Indeed, the decrease in {\color[rgb]{0,0,0}\sigma} results in a faster transition to the linear regime. The effect is particularly strong for {\color[rgb]{0,0,0}\sigma}=0.1, where most of the flow paths open very quickly after the first one. As we have discussed above, since the environment is more homogeneous, most of the possible paths are roughly equivalent.

4.4 Geometrical propeties of the flow field

We now study the statistical properties of the flow field. If the main quantity of interest should be the flow paths, their characteristics (branch length, rugosity, etc.) are quite difficult to define and measure precisely. Hence, it is more convenient to focus on the fluid at rest. We study the statistics of the non-flowing fluid clusters, i.e. the areas of fluid at rest that are surrounded by flowing channels. Using a standard Hoshen-Kopelman algorithm on the velocity map allows us to determine the cluster of fluid at rest. For each cluster, we then extract its size $\tilde{S}$. Its length {\color[rgb]{0,0,0}\tilde{L}} and width {\color[rgb]{0,0,0}\tilde{W}} are then determined by fitting it into a rectangle. Where the length, {\color[rgb]{0,0,0}\tilde{L}}, is defined as the dimension along the flow direction and the width, {\color[rgb]{0,0,0}\tilde{W}}, as the dimension transverse to it.

4.4.1 Probability distribution

We first study the size distribution. In Fig. 7 (left), we plot the cluster size distribution $P(\tilde{S})$ for a given $\sigma=1$ at different applied pressure difference. We can see that for any pressure difference, the distribution follows a power law for small sizes but with a cut-off at larger sizes. The value of this cut-off size decreases with the pressure difference. We then fit each distribution with the following law:

[TABLE]

where {\color[rgb]{0,0,0}\tilde{S_{0}}} is the characteristic cut-off size. As it can be seen in the inset of Fig. 7, the cut-off size is found to follow a power law with the flow rate:

[TABLE]

In Fig. 7 (right), we plot the same distribution but using rescaled variables according to the eqs. (20) and (21). We can then see that all the data are collapsing on a master curve. In addition, in Fig. 8, we traced the rescaled distribution for different heterogeneity parameter $\sigma$. As can be seen, the size of the clusters follows the same distribution function with the same exponents: {\color[rgb]{0,0,0}\tau} and {\color[rgb]{0,0,0}\gamma} for any sigma. We note, however, that the prefactor of the cut-off function varies with the disorder: {\color[rgb]{0,0,0}\tilde{S_{0}}}=D(\sigma)\;\tilde{q}^{-{\color[rgb]{0,0,0}\gamma}}, where D({\color[rgb]{0,0,0}\sigma})=3\sigma^{2}.

Eqs. (20) and (21) are thus reminiscent of what has been observed at the microscale. The power law distribution characterizes the multi-scale nature of the problem. At any pressure, cluster sizes are distributed within large range of sizes. The higher bound, {\color[rgb]{0,0,0}\tilde{S_{0}}} decreases with the flow rate (or pressure). Thus the closer we get to the global critical pressure, the wider the range of the power law distribution. Qualitatively, this decrease can be understood by the fact that the non-flowing areas are divided into smaller ones with the appearance of new branches as the pressure increases.

As for the flow rate - pressure cure, it is worth mentioning that if the general trend is similar to the miscroscale, the main difference lies in the exponent value. Indeed, for $\sigma=1$, we measure here {\color[rgb]{0,0,0}\tau}=1.15\pm 0.05 and {\color[rgb]{0,0,0}\gamma}=1\pm 0.05. The {\color[rgb]{0,0,0}\gamma} is thus similar the microscale ($\gamma=$ in [18]) but the exponent {\color[rgb]{0,0,0}\tau} is significantly different ({\color[rgb]{0,0,0}\tau}=1.5 at microscale). This indicates again that the process of branching is physically different for the two scales because the equations which are different (Stokes vs. Darcy) but also because the disorder is different (correlated log-normal field vs. random cylinder packing).

4.4.2 Cluster shape and Roughness

Now that we have shown that the cluster size is a multi-scale problem, we investigate the geometrical aspect of each cluster, as function of the applied pressure but also as function of the scale considered. For a given length scale {\color[rgb]{0,0,0}\tilde{L}}, we average the maximum width over clusters sharing the same length, \langle{\color[rgb]{0,0,0}\tilde{W}}\rangle_{{\color[rgb]{0,0,0}\tilde{L}}}. Moreover, we also study the shape of each cluster, $\tilde{w}(x)$, defined as the local width as function of the relative $x$ coordinate. We also average this width function over many clusters sharing the same total length: <\tilde{w}(x)>_{{\color[rgb]{0,0,0}\tilde{L}}}.

In Figure 9, we plot the contour shape \langle\tilde{w}(x)\rangle_{{\color[rgb]{0,0,0}\tilde{L}}} for different total lengths {\color[rgb]{0,0,0}\tilde{L}}. As expected, the cluster shape has a non-monotonic shape starting and ending with zero (e.g. \langle\tilde{w}(x=0)\rangle_{{\color[rgb]{0,0,0}\tilde{L}}}=\langle\tilde{w}(x={\color[rgb]{0,0,0}\tilde{L}})\rangle_{{\color[rgb]{0,0,0}\tilde{L}}}). The most interesting feature about the cluster shape is the fact that it is invariant to the scale as shown in Fig. 9. In this figure, we renormalize the contour shape by the maximum width and the $x$ coordinate by the total length {\color[rgb]{0,0,0}\tilde{L}}. All the contours then collapse on a parabolic-like master curve. This curve can be reasonably fitted to:

[TABLE]

with {\color[rgb]{0,0,0}\zeta}=0.75\pm 0.05. We also observed that the renormalized clusters’ shape (not shown) is also invariant with {\color[rgb]{0,0,0}\sigma} and can be fitted with the same exponent.

In addition, in Fig. 10, we plot the maximum width according to the length of the cluster for different {\color[rgb]{0,0,0}\sigma}. We observe that the maximum width also follows a power law:

[TABLE]

with an exponent independent of the disorder. The curve is, however, shifted with {\color[rgb]{0,0,0}\sigma}, reflecting the fact that the aspect ratio varies with {\color[rgb]{0,0,0}\sigma}. At low {\color[rgb]{0,0,0}\sigma}, the clusters are more elongated in the direction of flow, as can be clearly seen in the snapshot of the figure 3.

These results demonstrate the fractal (self-affine) nature of the flow field. At a given pressure, the size of the clusters is very widely distributed but with a shape which is invariant with the scale considered.

It is important to note that this behavior results from the self-affine roughness of the flowing paths surrounding the cluster. For instance, if one would assume that a cluster is delimited by two paths from a random walk model, as proposed by [38] in another context, the cluster shape would then be described by a first-time return random walk. It would thus follow relationships like eqs. (22) and (23) but with an exponent {\color[rgb]{0,0,0}\zeta}=1/2, which is the roughness of a random walk path.

Similarly as what has been observed earlier, the trend is comparable to the microscale. But here also, the main difference lies in the exponent value ({\color[rgb]{0,0,0}\zeta}=2/3\pm 0.03 at the pore scale) which is slightly higher. At the microscopic scale, the exponent could be justified with the roughness of the directed polymer problem.

We can put forward the following arguments to explain this small discrepancy. First, due to numerical constraints, our statistical average is relatively low (ten realizations). Secondly, the difference could also come from the log-normal distribution of the threshold, which might be considered as sufficiently ”extreme” to deviate from the standard directed percolation.

It can also be argued that the paths are not necessarily directed in our problem. But, as it can be seen in the snapshots of Fig. 3, for low {\color[rgb]{0,0,0}\sigma}, overhangs are no longer present but the roughness still remains the same.

5 Discussion/Conclusion

We have studied the flow of Bingham fluid in heterogeneous macroscopic porous media. We have shown that at the macroscppic level (Darcy’s scale) the flow properties share common features with the problem at the pore level. First, due to the disorder, the mean flow rate - pressure curve evolves non-linearly over a certain pressure range. Secondly, in this range, the velocity field exhibits multi-scale (fractal) properties that are reminiscent to other statistical physical problems with a critical transition (e.g percolation, avalanches, etc.). These multi-scale properties are described by power-law which are characterized by their exponent, {\color[rgb]{0,0,0}\alpha}, {\color[rgb]{0,0,0}\tau}, {\color[rgb]{0,0,0}\gamma} and {\color[rgb]{0,0,0}\zeta}. Yet, a significant difference between the two scales lies in the value of these exponents. At the macroscopic level, the mean flow rate increases with a power {\color[rgb]{0,0,0}\alpha} close to $3$ whereas at the microscale this exponent is close to $2$. We attribute this difference to the fact that at the macroscopic level, each individual flow path widen as the pressure increases. The exponents describing the flow geometry are also relatively different. At the macroscopic level, we observed ({\color[rgb]{0,0,0}\tau},{\color[rgb]{0,0,0}\gamma},{\color[rgb]{0,0,0}\zeta})=(1.15\pm 0.05,1\pm 0.05,0.75\pm 0.03), whereas at the microscopic level we had ({\color[rgb]{0,0,0}\tau},{\color[rgb]{0,0,0}\gamma},{\color[rgb]{0,0,0}\zeta})=(1.5\pm 0.05,1.1\pm 0.05,0.69\pm 0.03). Thus, the distribution of sizes ({\color[rgb]{0,0,0}\tau}) is very different, the cluster roughness is slightly different and the evolution of the cut-off size with the flow rate is similar. Another interesting results is the evolution of these scaling laws with the parameter of disorder {\color[rgb]{0,0,0}\sigma}. Indeed, despite the fact that the flow field is at first sight very different when we vary {\color[rgb]{0,0,0}\sigma}, we observed that the exponents remain the same. Only the prefactor of the different scaling laws are varying with {\color[rgb]{0,0,0}\sigma}.

We plan to extend this work in several directions. First, it would be interesting to investigate different types of permeability disorder. Even though the log-normal distribution is a very popular model in geostatitics it would be interesting to investigate other distributions. Since the branching mechanism results from the balance between finding the most favorable regions and the length to reach them, it should depend on the statistics of the large permeability events. One could thus expect to observe different scaling laws if the distribution presents fat tails (e.g power law). Another direction could be to investigate the transport phenomena (temperature, chemical species, etc.) associated with this type of flow. Since the flow structure is very complex, one would expect a non-linear transport properties depending on the applied pressure.

Acknowledgements

This work benefits from a PhD funding (R. Kostenko) from the French national radioactive waste management agency (Andra)