Does Rheology of Bingham Fluid Influence Upscaling of Flow through Tight Porous Media?

: Non ‐ Newtonian fluids may cause nonlinear seepage even for a single ‐ phase flow. Through digital rock technologies, the upscaling of this non ‐ Darcy flow can be studied; however, the requirements for scanning resolution and sample size need to be clarified very carefully. This work focuses on Bingham fluid flow in tight porous media by a pore ‐ scale simulation on CT ‐ scanned microstructures of tight sandstones. A bi ‐ viscous model is used to depict the Bingham fluid. The results show that when the Bingham fluid flows through a rock sample, the flowrate increases at a parabolic rate when the pressure gradient is small and then increases linearly with the pressure gradient. As a result, an effective permeability and a start ‐ up pressure gradient can be used to char ‐ acterize this flow behavior. By conducting flow simulations at varying sample sizes, we obtain the representative element volume (REV) for effective permeability and start ‐ up pressure gradient. It is found that the REV size for the effective permeability is almost the same as that for the absolute permeability of Newtonian fluid. The interesting result is that the REV size for the start ‐ up pressure gradient is much smaller than that for the effective permeability. The results imply that the sample size, which is large enough to reach the REV size for Newtonian fluids, can be used to investigate the Bingham fluids flow through porous media as well.


Introduction
Flow in porous media is ubiquitous in nature and industrial application [1][2][3].Darcy's law is the basic formula to depict the single-phase flow in porous media.As the Darcy's law can be derived from the Stokes equation, the Darcy's law is limited to low-Re Newtonian fluid flows in non-deformable porous media [4][5][6][7].Any deviation from the above conditions leads to non-Darcy flows [8][9][10][11][12], which have commanded great interest in the past twenty years.
The Non-Newtonian fluid [13,14] is one of the common reasons that cause a non-Darcy flow.Experiments have shown that waxy and heavy crude oils can be one of power-law fluids or Bingham fluids [15,16], which are two typical non-Newtonian fluids.The Bingham fluid through porous media can lead to a "start-up pressure gradient" [17][18][19][20] caused by the yield stress.The "real" start-up pressure gradient defines the pressure gradient below which there is no flow at all while the "pseudo" start-up pressure gradient is measured by fitting the data from the high-velocity flow regime.Experiments and simulations have been conducted to study the relation between flowrate and pressure gradient when a Bingham fluid flows through porous media.George et al. [21] and Bauer et al. [20] did measure the flow of yield-stress fluids through packed beads and found the power-law relation at the high-velocity flow regime.Talon and Bauer [22] performed numerical simulations of Bingham fluids through stochastically reconstructed porous media Citation: Liu, T.; Zhang, S.; Wang,

M. Does Rheology of Bingham Fluid
Influence Upscaling of Flow through Tight Porous Media?Energies 2021, 14, 680.https://doi.org/10.3390/en14030680 Academic Editor: Jens Birkholzer and found a parabolic transition flow regime and the linear flow regime when the pressure gradient increased.Chevalier et al. [23] also performed lattice Boltzmann simulations on 2D statistically reconstructed structures to study the statistical properties of Bingham fluid flows.
Because of strong backgrounds and important applications in oil/gas industry, the mechanism of Bingham fluid flows through real porous media has been studied at porescale by a powerful tool of digital rock analysis (DRA) [24][25][26][27][28]. DRA can investigate singlephase or multi-phase flows in real porous microstructures by simulating transport processes through reproduced pores from images by X-ray computed tomography (CT) or other imaging techniques.Bauer et al. [20] simulated a yield-stress fluid through sandstone and obtained the "true" or "pseudo" start-up pressure gradient.Yet, the disappointing fact is that the simulation results usually did not match the corresponding experimental data.On one hand, of course, the measurement technology needs improvements, especially for the low flowrate and long-term measurements.On the other hand, the accuracy of DRA needs to be seriously reconsidered and carefully treated.
Despite more and more advances in imaging techniques, there is still a contradiction between the imaging resolution and the sample size.A higher resolution means a higher possibility to recognize finer pores and pathways, which are very important for analysis of tight rocks.On the other hand, a higher resolution, at the current stage of imaging technologies, leads to a smaller sample size and a higher cost [29][30][31].For examples, nano-CT with a resolution of tens of nanometers requires the sample size to be tens or hundreds of microns.Therefore, we must find a balance or trade-off between scanning resolution and sample size.Although the higher-resolution imaging can output more details of real rocks, the corresponding smaller sample size risks losing representativeness and the results may be meaningless.It requires a careful consideration of both scanning resolution and sample size to achieve meaningful analysis results with appropriate costs.Studies have been done on the effects of sample size or scanning resolution, however most of them considered only one factor or both but separately [32][33][34][35][36][37][38].In our very recent work, the accuracy of DRA, which depends on the cut-off resolution (COR), sample size and the image processing method, has been quantified for the first time for Newtonian fluid flow in high-and low-permeability rocks [39].To our best knowledge, the upscaling of DRA with REV and resolution analysis for non-Newtonian fluid flow in tight rocks has never been reported.
In this study, DRA is used to study the Bingham fluid flow through tight porous media.The parabolic and linear flow regimes are obtained numerically, and one upscaling formula is proposed for the Bingham flow through porous media.The REV and COR of a "pseudo" start-up pressure gradient will also be investigated, and the results are compared with the Newtonian fluid flow.

Structure & Oil Properties
In this work, tight sandstone samples have been scanned with micro-CT.After checking the connectivity and REV of porosity of these samples, one sample is chosen as the target sample.The voxel size of the scanned images is ∆ 0.28 μm and the resolution is the reciprocal of voxel size  1 ∆ ⁄ 3.57 pixel μm ⁄ .A cube of  1000 is chosen as the initial box of region of interest and the image size  represents the voxels in each side of the cubic box.The raw structure is reconstructed following a typical image processing procedure, pre-processing, image segmentation and post-processing [40,41].In the pre-processing part, the brightness and contrast of the raw images are adjusted to make the images clearer, and the median filter is used to smooth the images and reduce the background noises.After that, a segmentation is completed with the Huang's algorithm [42] and the images are segmented into pores and solids.Finally, the connected pores are extracted in the post-processing part, which removes the tiny, isolated pores that are most probably produced during image processing and prepares for the flow simulations in the next step.
The previous study on this rock sample guaranteed that this sample satisfied the requirements of the representative element volume (REV) size and cut-off resolution (COR) for the permeability estimation of Newtonian fluid [39].To study the effects of sample size and resolution on Bingham flow, varying size and resolution of structures are numerically reconstructed from the raw cubic structure of  1000 (shown in Figure 1a).The sample size is changed with cropping operation by keeping the cubic center fixed and extracting a smaller cubic zone.Resolution is changed with coarsening operation used in Ref. [39].Resolution ratio of new structure is defined as    ⁄ Δ Δ ⁄ .Structures of varying sample size and structure resolution are shown in Figures 2 and 3.
(a) (b) Note that there are two different ways to choose sub-size samples when studying the REV effect.One is changing the sample size with the sample center fixed [39], which is adopted in this work.The other one is randomly choosing the sample positions [35].Advantage of the second way is that several sub-size samples can be chosen, and the coefficient of variance can be obtained to display the statistical deviation.However, the number of effective sub-size samples decreases with the increasing sample size, since the overlap becomes larger between samples.When the sample size does not reach REV, the sample number and the coefficient variance may provide some information of the mean property at that sample size, but this information seems meaningless because the sample is not representative.In our understanding, it seems unnecessary to choose several sub-size samples and calculate the coefficient of variance when investigating REV size from one raw sample.The sample-dependent or position-dependent effect can be studied with several raw samples and REV needs to be satisfied for each sample.Thus, the first way is preferred to study the REV size from one raw sample and this way also matches the REV definition proposed by Bear [1].Experiments on fluid properties have shown a non-linear relationship between shear stress  and shear rate  [15] for non-Newtonian fluids.It is found that Bingham's model can be used to depict this non-Newtonian waxy crude oil property [14].In the present study, a bi-viscous model is used, which gives: to depict a linear relationship between  and  with an interception  . is the zeroshear rate viscosity and  is the infinite shear rate viscosity, which is much smaller than  . is the critical shear rate above which the effective viscosity decreases from  to  when  increases.In the model fit, the experimental data gives the model parameters listed in Table 1.For a Bingham fluid, it behaves like solid when it is not moving or moves slowly with very large viscosity but can also flow like a fluid with small viscosity when velocity gets larger.Toothpaste and some crude oil are typical examples of Bingham fluids.Though the fluid is to move theoretically even when  is not reached using this bi-viscous model, it is still a good representative of the Bingham fluid, especially when studying "pseudo" start-up pressure gradient in the high-flowrate regime, as shown in Figure 5.

Pore-Scale Modeling of Bingham Fluid Flow
The non-Newtonian fluid flow through porous media is simulated using the finite volume method (FVM) to solve the Navier-Stokes (NS) equations at pore-scale.The computation is performed on OpenFOAM [43], an open source c++ library for computational fluid dynamics (CFD).To ensure grid independence, the extrapolation flowrate calculation approach is taken from the numerical results at the reference and refined mesh [39].For all cases, the structured cubic lattices are adopted.Since the rock structures are obtained by CT images, we reproduce the structures by replacing every pore voxel with a cubic lattice.All simulations are terminated when the flow rate gets stable and the results are output at the steady state.
In our previous study [39], the Newtonian fluid flow through body-centered cubic (BCC) packing of spheres was simulated, which verified the numerical solver for Newtonian fluid flow in porous media.To examine the solver capability for the Bingham fluid flow, a Bingham fluid flow through 2D slit is considered.The simulation domain size is 200 100 as shown in Figure 6a.The fixed-wall condition is applied at upper and lower boundaries and pressures are given at inlet and outlet.Parameters for the bi-viscous model are given in Table 1.Theoretical analysis gives the velocity profile as: where  ⁄ gives the location that the effective viscosity changes to  .Figure 6b shows the velocity profile by numerical simulations, which agrees well with the analytical solution.The integral of velocity will give the formula of flow rate as: This means that the flow rate changes from a parabolic to a linear form when the pressure gradient increases.It also gives the "real" and "pseudo" start-up pressure gradient as  ℎ ⁄ and 3 2ℎ ⁄ .For 3D cases of flow in cylinder, the theoretical analysis on flow rate gives similar results as:

Upscaling of Bingham Fluid Flow
After validating the numerical methods, we consider the non-Newtonian fluid flow through porous media.Three-dimensional porous structures are reconstructed as shown in Figure 7. Flow is simulated through the pores.A no-slip boundary condition is adopted on the solid walls.Pressures are given at inlet and outlet so that a steady-state flow rate can be calculated.By varying the pressure difference between inlet and outlet, we can obtain the relationship between flowrate and pressure gradient.The Reynolds number is checked firstly at the largest pressure difference and it is found that    ⁄ ~10 , which means that the inertial effect makes little significance in this flow regime [44].Since the relation changes from parabolic to linear form, by assuming a continuous transform between two forms, the following formula is proposed: Here,  and  are the length and the cross-section area of the sample, respectively,  is the effective permeability and  is the "pseudo" start-up pressure gradient.This form differs from the relation of above 2D or 3D channel since the "real" start-up pressure gradient for this sample is small compared to the "pseudo" start-up pressure gradient.To estimate the precise "real" start-up pressure gradient demands more data points in the very-low pressure gradient regime, but we focus more on the large pressure gradient regime and the "pseudo" start-up pressure gradient in this work.For simplification, we omit "pseudo", and the following start-up pressure gradient means the "pseudo" start-up pressure gradient.This form satisfies the parabolic form between  and Δ  ⁄ when Δ  ⁄ is small and the linear form when Δ  ⁄ is large.It also guarantees a continuous change between these two forms.The results in Figure 8 shows that this model can well match the simulation results.One thing to mention is that this form does not fit the case where "real" start-up pressure gradient is not negligible when Bingham fluid flows through the sample.In that cases, one needs to change   ⁄ with   ⁄  , which means the distance to "real" start-up pressure gradient.More discussion and detail derivations on this upscaling formula can be found in the Appendix A. Because the mesh size may introduce numerical error to simulation results, the extrapolation scheme is used to exclude this error.As shown in Figure 9, the results change when the mesh size changes.However, the mobility (or flow rate) changes linearly with mesh size at varying pressure difference and the extrapolation result can exclude this numerical error and obtain more accurate results [39].With an extrapolated flow rate at varying pressure difference,  and  at varying structure size and the resolution can be accurately obtained.

Results
The Bingham fluid flow through porous structure is simulated for different pressure difference and different sample size and scanning resolution.The  and  are obtained from the relationship between  and Δ.The effects from fluid properties on the REV size and critical scanning resolution will be therefore studied.

REV Size
The relation between flow rate and pressure difference at varying structure size is shown in Figure 10.At least four high pressure-drop results are presented to guarantee accurate estimation of the effective permeability and the start-up pressure gradient.
Figure 11 shows that  and  changes with the sample size.For the effective permeability, a larger deviation is observed at smaller structure size.The effective permeability becomes stable when the structure size gets larger than 800, which is approximated as the REV size for  .As for the start-up pressure gradient, through deviation can be observed, the deviation is quite small compared to its mean value and a REV size of 100 is obtained for the start-up pressure gradient.This means that the start-up pressure gradient is a much more stable property compared to the effective permeability.

Cut-Off Resolution (COR)
The relation between the flow rate and the pressure difference at varying resolution is shown in Figure 12.At least four high pressure results are simulated to guarantee accurate estimation of the effective permeability and the pressure gradient.
Figure 13 gives the results that  or  changes with resolution.For the effective permeability, large deviation is observed at low resolution and it gets stable when the resolution gets higher than 0.375, which is defined as COR for  .As for the start-up pressure gradient, a similar trend is observed and the COR for start-up pressure gradient is also around 0.375.

Discussion
The REV size for Bingham fluid flow in porous media is an important issue for both theoretical analysis and industry applications.In some previous studies, it was usually believed that the REV size for non-Newtonian fluids is larger than that for Newtonian fluids.Theoretically, no significant difference has been found to estimate REV of effective permeability for Bingham-fluid and Newtonian-fluid flows.Our framework provides a tool to quantify REV and permeability for different types of fluid in porous media.Figure 14 shows the comparison of permeability estimated for Bingham and Newtonian fluid flows at varying structure size, which indicates that the Bingham fluid appears similar as the Newtonian fluid in permeability evolution.In high-velocity flow regime, this result seems trivial since the influence of yield stress on the effective viscosity is negligible.Meanwhile, one more parameter  needs to be considered for the Bingham fluid.Typical expectation of REV for  is larger than REV for permeability and thus the REV for Bingham fluid flow would be larger than that for a Newtonian fluid.However, the numerical simulations have given an opposite result that the REV for start-up pressure gradient is smaller than REV for  .
The theoretical analysis of Bingham fluid flow through cylinder indicates that  is proportional to  ℎ ⁄ while  is related to ℎ .This means that  is sensitive to changes of reciprocal of pore size and  is sensitive to changes of square of pore size.Resolution R The simulation results imply that the REV for 1 ℎ ⁄ is smaller than that for ℎ .Still, the reciprocal of specific area is used to indicate the average pore size, ℎ 1 , changes of 1 ℎ ⁄ and ℎ with structure size  can be plotted in Figure 15.The results give the estimate of the REV for 1 ℎ ⁄ and ℎ as 200 and 400 respectively.This confirms the above implication and explains why the REV for  is smaller than that for  .Therefore, contrary to the typical expectation, the REV for a Bingham-fluid flow is almost the same as that for a Newtonian-fluid flow.This means that the sample size, which is large enough to estimate permeability of Newtonian fluids, satisfies the REV requirement to study Bingham-fluid flow through porous media.

Conclusions
Digital rock analysis provides a powerful tool to investigate the transport of complex fluids in porous media while the critical size of REV and scanning resolution should be satisfied to guarantee representativeness and accuracy.In ] Structure size N resolution, the numerical framework obtains the REV and COR for  and  .It is found that the COR of  is similar with that of  but the REV for  is much smaller than that for  .The present study indicates that the REV and COR for Bingham fluid is almost the same with those for Newtonian fluid.Therefore, the samples, which satisfy REV and COR requirements for Newtonian fluid, can also be used to study the Bingham fluid flow and to estimate the start-up pressure gradient in tight porous media.
When  approaches , the flow rate is  *  1 .

(A7)
For large viscosity ratio,  approaches 1 and  * approaches 0. Define   , the flow rate near  is written as , flow rate is re-written as: For large viscosity ratio,  approaches 1 and  With  * 1, pressure gradient needs to satisfy condition  2.For Bingham fluid, this condition reflects the existence of start-up pressure gradient.
When  is large enough, flow rate can be approximated as: This is the linear relation form for large pressure gradient.When  approaches 2, the flow rate is:

(A17)
For large viscosity ratio,  approaches 1 and  * approaches 0. Define  2 , the flow rate near 2 is written as , flow rate is re-written as:

Figure 1 .
Figure 1.3D visualization of pore structure of Changqing tight sandstone sample and its pore size distribution.(a) Pores and solids are both shown.Pores are shown in black color and solids are shown in gray color.(b) Pore size distribution of tight sandstone with characteristic pore diameter  1.4 μm.

Figure 2 .
Figure 2. Varying size structures obtained from cropping operation.The resolution ratio is  1.0 and sample size changes from  1000 to  250 with a step of 250.Pores are shown in black and solid in white.Only a 2D slice is presented in the figure.

Figure 3 .
Figure 3. Varying resolution structures obtained from coarsening operation.The sample size is  1000 and resolution ratio changes from  1.0 to  0.25 with a step of 0.25.Pores are shown in black and solid in white.Only a 2D slice is presented in the figure.

Figure 4
Figure 4 presents the porosity and reciprocal of specific area for varying structure sizes.Reciprocal of specific area is used to represent the averaged pore size.It is found that REV of porosity and reciprocal of specific area reaches  200, which guarantees the representativeness of the rock structure in the following analysis of flow property.

Figure 4 .
Figure 4. Results of (a) porosity  and (b) reciprocal of specific area  change with structure size .

Figure 5 .
Figure 5. Experiments data and the model fit of the crude oil viscosity.(a) Experiments data of the crude oils: shear stress vs. shear rate relation [15].(b) Bingham model fit of one typical crude oil viscosity experiment data from (a).

Figure 6 .
Figure 6.Validation on Non-Newtonian fluid flow through a 2D slit.(a) Schematic of the 2D parallel plates.(b) Velocity profiles by numerical simulations and analytical solution.

Figure 7 .
Figure 7. Computational domain of reconstructed 3D porous structure at  1000 and  0.25.Only pore part structure is presented.The color reflects the pressure distribution.

Figure 8 .
Figure 8. Relationship between flow rate and pressure difference for Bingham fluid flow by simulation.The sample size is  100 and resolution ratio is  1.0.The proposed model (dashed line) matches the simulation results (dot) well.The flow rate changes from a parabolic form to a linear form when the pressure difference increases.

Figure 9 .
Figure 9. Numerical results for Bingham fluid flow simulation at varying mesh size and pressure gradient.The sample size is  100 and resolution ratio is  1.0.(a) Relationship between flow rate and pressure difference at varying mesh size.The index  in the figure represents the refinement level. 2 means that each voxel mesh is divided into 2 8 reined meshes.(b) Mobility (or flow rate) changes linearly with mesh size for varying pressure difference simulation.The index  in the figure represents the pressure difference applied in the simulation.

Figure 10 .Figure 11 .
Figure 10.Numerical results for Bingham fluid flow simulation at varying pressure gradient ∇ and structure size .The resolution ratio is  1.0.(a) Relationship between mean velocity and pressure gradient at varying structure size.(b) Same results shown in (a), but pressure gradient range changed to keep only 4 data points in the linear form regime.The index  in the figure represents the structure size.

Figure 12 .Figure 13 .
Figure 12.Numerical results for Bingham fluid flow simulation at varying pressure gradient ∇ and resolution .The sample size is  1000 and resolution ratio changes from  1.0 to  0.125.The index  in the figure represents the resolution ratio.

Figure 14 .Figure 15 .
Figure 14.Results of permeability change with structure size for Bingham and Newtonian fluid flow.
this work, these critical requirements for non-Newtonian fluid flow through tight porous media are examined.Porous structures with various size and resolution are reproduced to study REV and COR for single-phase Bingham fluid flow and the results are compared with the ones for Newtonian fluid flow.The results show that the relation between flow rate and pressure gradient for Bingham fluid flow changes from a parabolic form to a linear one as the pressure gradient increases.An upscaling model is proposed connect the effective permeability and the start-up pressure gradient to depict this flow.Considering varying sample size and 0

)
This is the parabolic relation form.To show the validity of flow rate model proposed in Section 2.3, the results of analytical solution and model fit in shown in FigureA1.It can be seen from the figure that the model can match the analytical solution well.

Figure A1 . 1 .
Figure A1.Results of non-dimensional flow rate changes with non-dimensional pressure gradient for Bingham flow through 2D parallel plate.The flow rate changes parabolically when pressure gradient is small and linearly when pressure gradient is large.The model that links parabolic and linear form can well match the analytical solution obtained.Appendix A.2. Bingham Flow in 3D CylinderAnalytical solution gives the velocity distribution as     ℎ  ℎ , 0    ℎ  ℎ ,   ℎ parabolic relation form.To show the validity of flow rate model proposed in Section 2.3, the results of analytical solution and model fit in shown in Figure A2.It can be seen from the figure that the model can match the analytical solution well.

Figure A2 .
Figure A2.Results of non-dimensional flow rate changes with non-dimensional pressure gradient for Bingham flow through 3D cylinder.The flow rate changes parabolically when pressure gradient is small and linearly when pressure gradient is large.The model that links parabolic and linear form can well match the analytical solution obtained.

Table 1 .
Bingham fluid model parameters for an experimental crude oil.