Comparative study of high-pressure fluid flow in densely packed granules using a 3D CFD model in a continuous medium and a simplified 2D DEM-CFD approach

An isothermal compressible single-phase fluid flow through a non-homogeneous granular body composed of densely packed overlapping spheres imitating rock under high pressure was numerically studied using two different approaches. The first approach called the full 3D CFD model used the finite volume method (FVM) to solve the Reynolds-averaged Navier–Stokes equations using Reynolds stress model (BSL) in the continuous domain between the granulates. The model was verified, based on experimental and numerical results from the literature. The second approach was a simplified coupled DEM-CFD model based on a fluid flow network. The main aim of the work was to develop a validation procedure for simplified coupled DEM-CFD models due to the lack of experimental data for fluid flow characteristics in densely packed granules under extremely high-pressure conditions. First, a series of numerical simulations were performed for the fluid domain with the full 3D CFD model. The results of those simulations were next used to validate the 2D numerical results of the simplified coupled DEM-CFD model with respect to velocities, pressures, densities and flow rates. Almost the same pressure and density distributions and mass flow rates were obtained in both approaches. However, the fluid velocity was different due to the different fluid volumes in both fluid domains. The current simulation results constitute a reliable benchmark for validating other coupled 2D/3D DEM-CFD models that use a fluid flow network approach.


Introduction
There are several hydro-mechanical engineering problems [1][2][3][4][5][6][7][8] in which it is necessary to study physical phenomena using the discrete element method (DEM) coupled with computational fluid dynamics (CFD). This coupling is of major importance for analyzing mesoscale fluid-solid interactions in porous heterogeneous materials (e.g. rock and concrete with porosities significantly below 10%) that strongly affect their global behaviour (particularly when different failure patterns appear). As compared with usual continuum methodologies, the coupled DEM-CFD approaches are more realistic since they allow for a direct simulation of the material meso-structure and thus are very useful for studies of mechanisms of the initiation, growth and formation of cracks and fractures.
Recently, several simplified coupled approaches combining the discrete element method (DEM) and computational fluid dynamics (CFD) (based on fluid flow networks) were formulated mainly in 2D conditions for describing hydromechanical problems (e.g. hydraulic fracturing process in rocks [9][10][11][12][13][14][15][16][17][18][19]). They consider incompressible/compressible laminar viscous fluid flow and interaction mechanisms between flowing fluid and particles with the so-called pore networks, called the pore-scale finite volumes that are built through a weighted Delaunay triangulation over the discrete element packing. The edges of triangles connect the gravity centers of discrete elements. The different numerical methods are applied to solve the governing equations of motion. The most common is a single-phase fluid flow model in pores and cracks/capillaries (full saturation); the only phase is a liquid. If numerical investigations consider the material behaviour with e.g. rock-like properties, other numerical methods such as FEM, FVM, LBM and SPH cannot be easily coupled to DEM that requires densely packed elements. This requires then extremely dense meshes due to both enormous pressures and velocities in fluid or a complex fluid domain geometry (e.g. during a hydraulic fracturing process in rocks). Meso-structures of some engineering materials (e.g. rocks or concretes) possess a very low porosity that leads to a very high packing of discrete elements. Moreover, meso-scale processes are often unsteady (time-dependent) and characterized by pronounced changes in the geometry topology what contributes to frequent remeshing processes (even computing parallelization cannot solve this problem [20]). The very densely packed discrete domains cause the following issues for classical CFD models: An additional difficulty in using classic CFD models to study mesoscale phenomena in materials with very low porosity (10% and less) is the lack of experimental results (e.g. during hydraulic fracturing). In the practice, DEM-CFD models can only be calibrated for the bulk permeability that is a material property to be experimentally determined. However, this does not guarantee that the fluid flow is realistically reproduced at the meso-level. Considering that the purpose of meso-scale investigations is to understand local physical phenomena affecting the material behavior at the macro-scale, this validation method may also be insufficient. Hence, the only method of validating the simplified DEM-CFD model and the characteristics of fluid flow in the material at the mesoscale is a comparison of calculation results with corresponding outcomes for a continuous domain between discrete elements (not reduced to a system of channels) with commonly used CFD models. While the overall trend of pressure distribution in the media is expected to be almost linear, an incorrect DEM-CFD model may not even be able to reproduce this trend. Therefore, the validation of DEM-CFD models should be performed on very simple specimens with known macroscopic and mesoscopic fluid flow properties. The best practice is to prepare a specimen with a defined porosity with boundary conditions reproducing a one-dimensional fluid flow at the macro-scale. The specimen porosity must be as low as possible but still allowing for mesh generating.
The main research goal is to propose a validation procedure for simplified 2D fluid flow models coupled with DEM. The validation procedure is presented with our coupled 2D DEM-CFD approach that was formulated and used for describing a hydraulic fracturing process in rocks [21,22]. It simplifies the flow of fluid in a continuous domain under isothermal and high-pressure conditions to laminar flow through a virtual pore network composed of channels between granulates (VPN). It considers the solid domain (discrete elements) and the fluid domain at the same time. It was significantly improved as compared to the existing fluid flow network models [9][10][11][12][13][14][15][16][17][18][19]. The approach improvements focused on the detailed tracking of liquid and gas fractions in pores and fractures depending on their different geometry, size, and position under isothermal conditions [21]. In addition, our approach was extended to multiphase fluid flow [22]. The coupled DEM-CFD approach was developed by the authors and implemented into the opensource DEM program YADE [23,24]. The approach was directly compared with numerical simulation results of fluid flow in a continuous 3D domain between very densely packed and overlapping spheres (imitating a rock specimen) under high pressures (ca. 70 MPa) on isothermal conditions, based on the Reynolds Averaged Navier-Stokes equations using the Reynolds Stress Baseline (BSL) turbulence model (hereinafter referred to as the full CFD 3D model). In this model, the fluid domain between discrete elements was considered only. The full CFD simulations of fluid flow were carried out with the commercial software ANSYS CFX [25]. The full CFD approach was validated with an experiment [26] and its numerical simulation [27]. The velocities, pressures, mass flow rates, and densities in a non-homogeneous densely packed domain with the same geometry and permeability were directly compared for both the approaches able to simulate fluid flow. No fractures were taken into account in a deforming non-homogeneous granular specimen (they will be considered in the next research step). The current simulation results constitute a reliable benchmark for validating other coupled 2D/3D DEM-CFD models that use a fluid flow network composed of channels in a domain between spheres (such a benchmark calculation has not yet been demonstrated in the literature). Another goal of our present investigations is to compare the calculation time in 2D and 3D simulations to present a numerical model with a feasible computation effort.
The paper is arranged as follows. After Introduction (Sect. 1), the fluid flow models using the full 3D CFD and simplified 2D DEM/CFD are described in Sects. 2 and 3. The numerical validation of the fluid flow model using the full 3D CFD in a continuous domain between spheres, based on an experiment and its numerical simulation, is discussed in Sect. 4. Section 5 reports on comparative numerical study results of the fluid flow in a densely packed granular specimen imitating rock using both the full CFD and simplified DEM-CFD approach. Finally, some concluding remarks are summarized in Sect. 6.

Full CFD model of fluid flow in continuous 3D domain between spheres
Turbulent compressible single-phase fluid flow through a heterogeneous very dense packing of fixed granules (imitating a rock specimen) was numerically studied under isothermal conditions. Although spherical particles with the same diameter were used, the final shape of each particle was not spherical due to overlaps between them. It was also assumed that there was no mass source or external body force in the fluid domain. Gravitational forces were neglected. An unsteady turbulent model was assumed. The results were presented at the equilibrium state.
In fluid mechanics, governing equations are derived from conservation laws: conservation of mass, momentum, and energy. The conservation of energy equation is not considered in this study due to isothermal conditions. The mass conservation equation known also as the continuity equation is given by where ρ is the fluid density, t is the time and V is the velocity. The first term presents the density variation in time in a control volume and the second term presents the mass flow rate from one surface to the other surface of the control volume. The density varied according to a barotropic model, defined as where 0 is the density for the reference pressure, P 0 is the reference pressure and C is the fluid compressibility for water at 70 • C ( 0 = 977.36 kg/m 3 , P 0 =0.1 MPa and C = 4e −10 1/Pa). The momentum conservation equation for compressible flow is where P is the fluid pressure and is the dynamic viscosity. Equation 3 describes the momentum transport under assumed conditions for both laminar and turbulent flow regimes. For distinguishing a dominant flow regime, the Reynolds number Re is used that defines the ratio of inertial forces to viscous forces where L is the characteristic length scale. However, when a fluid flows through a porous medium, Re is usually related to the average particle diameter and porosity [28]. If a porous structure is reproduced by an infinite number of square bars/ tubes, the Reynolds number Re pore is used based on the pore size [29]. In the case of DEM, the material structure is reproduced by discrete elements, e.g. spheres, and the following definition of the Reynolds number Re part is usually used [30] where D p is the average particle diameter, ϕ denotes the porosity ( = volumeoffluid totalvolume ) , is the dynamic fluid viscosity,Vy is the pore average streamwise velocity and U s denotes the superficial velocity (U s = Q/A t , where Q and A t are the volumetric flow rate and total bed cross-sectional area, respectively). The high-velocity flow in porous media leads to microscopic (pore level) turbulences [31,32]. Highly chaotic structures usually develop when Re part > 120 [30] or Re pore > 300 [32]. The mean diameter of spheres is used as a characteristic length in the particle Reynolds number and the mean pore size is used in the studies of porous media composed of tubes. The Re number can be also defined as the permeability Reynolds number Re k = Us √ K D , where K D denotes Darcy's permeability coefficient. This Re number, Re k , has no physical meaning and it should not be confused with the ratio between inertia and viscous forces (Eq. 4) [31]. The permeability coefficient K D is measured using Darcy's law that constitutes a linear relationship between the fluid velocity and pressure drop: where ∇P is the pressure drop over a given distance. Darcy's law is solely valid up to the creep flow limit with Re k < 0.1 [28,31] (the pressure drop within the medium is solely due to the microscopic viscous drag). As the flow velocity increases, a correction term known as the Forchheimer term is added to Darcy's law to consider the drag due to both inertial effects inside pores and turbulent dissipation [32]: where K F is the Forchheimer permeability coefficient that is close to but not the same as K D [33]. The flow in porous media has two main differences from the free-stream turbulent flow: the size of eddies is limited by pore size and additional drag is caused by a porous solid matrix [32]. The turbulent flow was analyzed in the current paper with the Baseline (BSL) Reynolds Stress model that solves two transport equations: one for the turbulent kinetic energy k and one for the turbulent frequency [25]. The BSL model is based on the ω-equation which precisely describes the fluid behaviour near walls. In BSL, an equation for the transport of all components of Reynolds stresses tensor is used without the eddy viscosity hypothesis [25]. The Reynolds Stress model directly introduces the 'Reynolds stresses' terms that are computationally expensive due to six additional equations. The BSL model is better suited to complex flows such as fluid flow in specimens used in this study. This is due to a precise model of the turbulence production terms ( P k , P kb and P b in Eqs.8 and 9) and the model of stress anisotropy [25]. where is the turbulent frequency, k denotes the turbulent kinetic energy (defined as 0. ′ are the velocity fluctuations in the OX, OY and OZ direction, P k is the production rate of turbulence (turbulence due to viscous forces), P kb and P b are the buoyancy production terms that represent the influence of buoyancy forces on kinetic energy and frequency, respectively and σ ω2 is the turbulent Prandtl number in the -region. U is the average velocity component, F 1 is the blending function equal to 1 near the surface; it decreases down to 0 outside the boundary layer, t is the eddy viscosity or turbulent viscosity (defined as t = k ) and ′ denotes the model con- .09). The turbulent Prandtl number k and and the coefficients and are mixed between values from two sets of constants, corresponding to the constants for the -based model and the constants for -based model transformed to the -formulation [25]. The constants are k1 = 2, 1 = 2, 1 = 0.075, 1 = 0.553 for the Reynolds stress-the -turbulence model in the -region and k2 = 1, 2 = 1∕0.856, 2 = 0.0828, 2 = 0.44 for the Reynolds stress-the -turbulence model in the -region [25].

Virtual pore network model of coupled DEM-CFD
The DEM-CFD model couples the discrete element method in a discrete domain with computational fluid dynamics in a continuous domain. DEM directly simulates the material meso-structure and thus is suitable for comprehensive studies of mechanisms of the initiation, growth and formation of localized zones, cracks and fractures at the meso-scale that greatly affect the macroscopic behaviour of frictional-cohesive materials [34,35]. It easily represents discontinuities caused by fracturing or fragmentation. The ease of access to information at the particle scale makes DEM a very useful tool to study the dynamics of systems composed of particles. The DEM disadvantage is a huge computational cost. The DEM calculations were performed with the three-dimensional spherical explicit discrete element open code YADE [23,24], allowing a small overlap to be formed between two contacted bodies (soft-particle model). An arbitrary micro-porosity can be achieved in DEM since the particles may overlap. The DEM model used is briefly summarized in Appendix 'A'. The general concept of a fluid flow algorithm using DEM was adopted from Cundall [36], Hazzard et al. [37], and Al-Busaidi et al. [38]. In this concept, fluid flow was simulated by assuming that each particle contact was an artificial flow channel (between two parallel plates in 2D or along a duct in 3D) and those artificial channels connected real reservoirs in particulate media (pores, fractures, and pre-existing cracks) that stored fluid pressures. Thus, the pressure in reservoirs depended both on the mass transported along channels from/ to other reservoirs and the volume changes of reservoirs. Since the volume of reservoirs changed due to the material deformation (described by discrete elements in DEM), the fluid density had to also change (the fluid in reservoirs must be compressible). Thus, the fluid moved in channels while the reservoirs solely stored pressure. The artificial channels created a fluid flow network. The fluid flow in artificial channels was characterized by a simplified laminar flow of the incompressible fluid instead of a compressible fluid model in real reservoirs. However, the spatially constant density of the fluid in channels might change in time due to density changes in reservoirs.
The original system consisted of two coexisting domains: 3D discrete domain (spheres) and 2D fluid domain. The gravity centers of spheres were located on the XOY specimen mid-plane. To create a two-dimensional fluid domain, spherical 3D particles were projected onto the mid-plane domain. Next, a remeshing procedure discretized the overlapping circles (projected spheres), determined the contact segments and deleted the overlapping areas [21]. All discretized pores between discrete elements called here virtual pores (VP), created a virtual pore network (VPN) to precisely reproduce their changing geometry (shape, area, and location). The displacements of spheres in the perpendicular direction OZ and the rotations around the axes OX and OY were fixed.
The vertices of triangles were located on circle circumferences. The gravity centres of grid triangles (VPs) were connected by channels composed of two parallel plates that created a fluid flow network. The isolated pores in 2D were not isolated in 3D. The virtual (artificial) channels were introduced in the 2D fluid flow network to reproduce real flow in 3D. They were located between spheres in contact and connected the isolated pores. Hence, two types of channels were introduced ( Fig. 1): (1) the channels between spheres in contact and connecting the isolated pores (the so-called virtual channels S2S) ( Fig. 1a) and (2) the channels connecting grid triangles in pores (T2T) (Fig. 1b). The channels 'T2T' and 'S2S' connected the gravity centres of triangles and created a fluid flow network (VPN) (Fig. 2). The hydraulic aperture h of virtual channels (S2S) was related to the normal stress by a slightly modified empirical formula of Hökmark et al. [40]. The hydraulic aperture of the channel type T2T was directly related to the geometry of adjacent triangles. In general, all particles could move during simulations due to the applied fluid pressure. Consequently, the contact stresses and the overlap areas changed. The fluid flow network was then modified. If particle bonds were not broken, the calculations were carried out with 'T2T' and 'S2S'. If bonds were broken, the virtual channels 'S2S' that connected two separated pores by overlapping particles were replaced by the channels 'T2T'. The new fluid flow network was generated by a remeshing procedure for actual particle positions. The geometric criterion was used to trigger the remeshing procedure that was switched on if the sum of the particle radii was less than the distance between the particle centers of gravity.
The fluid moved through channels through a thin layer separated by two parallel plates close to each other. Virtual pores accumulated pressure and stored density. The fluid density change in VPs was related to the mass change that resulted in pressure variations. The equation of momentum conservation was thus neglected in triangles but the mass was still conserved in the entire volume of triangles. The equations of state and continuity were employed to compute the density of fluid stored in VPs. The fluid flow in channels was estimated by solving continuity and momentum equations for the laminar flow of incompressible fluid.
For single-phase flow, the mass flow rate in the channels was estimated, based on the Poiseuille flow law [39], derived from the classical Reynolds lubrication theory [41]: where This algorithm was repeated for each VP in VPN using an explicit formulation (Fig. 3). The discretization algorithm was based on the alpha-shape theory [42] and the Delaunay triangulation. The grid remeshing was automatically performed when the topological properties of the grid geometry changed [21]. The computational results (e.g. pressures) were accurately transformed from the old grid to the new one, by assuming that the mass was a topological invariant. The coupling scheme of DEM with CFD (described in detail in [21]) involved two sets of discrete equations to be solved: the flow rule defined for all VPs and the law of motion in DEM for all discrete elements. The two-way coupling scheme was based on a transfer of pressure and shear stress forces from CFD to DEM and the time derivative of VP volumes from DEM to CFD [21]. The pressure and shear forces from the fluid caused sphere displacements in DEM that changed coordinates and volumes of triangles in VPN. The time derivative of the VP volume was computed in DEM and then transferred to CFD (Eq. 11). As a result, the volume change in DEM affected the pressure change in the fluid. The fluid pressures in VPs and 'S2S' channels were subsequently converted into the forces acting on spheres. The pressures in VPs and shear stresses in 'S2S' channels were calculated by CFD in each time step and transferred from CFD to DEM. They were used to compute fluid forces that were added to contact forces before the time integration to update the displacements of each discrete element. The shear stresses in 'T2T' channels were also converted into forces acting on spheres.
In the DEM-CFD approach, the density of the Delaunay triangulation grid depends on the number of segments into which the circle (sphere cross-section) is divided. Consequently, the number of channels in the fluid flow network depends on the number of divisions. The griddependency tests were performed. For a division number greater than or equal to 12, the results changed by less than 5%. Thus, a division number equal to 12 was assumed in all simulations.
The pure CFD was calibrated with the aid of permeability laboratory tests on rock specimens [22]. The following material constants were needed for the CFD simulations: reference pressure P 0 , fluid density 0 for the reference pressure P 0 , fluid bulk modulus K = 1 C , dynamic viscosity of the liquid μ. The inverse of the bulk modulus of the fluid gave a fluid compressibility C. From one layer of spheres used in a 2D problem, it is easy to extend the model into 3D calculations. All numerical parameters may be used without the necessity of re-calibration.

Experiment
Initially, to validate our full CFD model, the experimental results by Suekane et al. [26] on a packed bed of spheres were used for comparative purposes (spheres were fixed and arranged at each corner of the cube). A special apparatus was constructed to measure the velocities of the fluid propagating through voids of a packed specimen (Fig. 4). The velocities were measured with the magnetic resonance imaging (MRI) technique. The experimental apparatus included a flow channel aligned vertically inside the MRI system.
The diameter of the channel was d = 50 mm. The fluid channel was composed of 5 unit cells of the same size. Each unit cell 28 × 28 × 28 mm 3 had a cubic shape and contained eight fixed 1/8 spheres with a diameter D P = 28 mm. The spheres were placed in two layers with four 1/8 spheres. The volumetric porosity ϕ of the unit cell was ϕ = 47.6%. The circulating fluid was water supplied from a reservoir with a constant pressure to keep the flow rate constant. Three components of the fluid velocity vectors in the 3D void space were directly visualized. The y-component of the fluid velocity was registered in the fourth unit cell at y = 0 (Fig. 4) for five different Re part numbers from the inertial flow region to the unsteady laminar flow region. The velocities were presented along the x-axis (blue line in Fig. 4). The phase encoding method was used [43] to find the velocity distribution in a specimen cross-section. The Reynolds numbers Re part for the streamwise velocity for five arbitrarily chosen particles were: 12.17, 28.88, 59.78, 105.5, and 204.74. The dominant flow regime for Re part = 204.74 was reported to be laminar ( [27,[44][45][46]).
For the value of Re part = 12.17, the fluid velocity decreased with the increasing area of cross-section normal to the main flow direction (the local geometry significantly influenced the flow). The water flowed over the surfaces of spheres like a creeping flow in which viscous forces were dominant. When increasing Re part , the local geometry was not the only factor affecting the velocity variation due to an increasing influence of inertia forces over viscous forces. With the high value of Re part (e.g. 204.74), where the inertial forces were dominant, the fluid moved through the pores like a jet without a velocity change. The maximum velocity in the center of pores increased with a growth of Re part even four times as compared to the average velocity in pores. The velocity profile was always parabolic.

Numerical results
Our numerical results were also compared with similar numerical simulations of the experiment [26] performed by Gunjal et al. [27] with the commercial software ANSYS Fluent. Unstructured tetrahedral grids were assumed in their study to simulate incompressible laminar flow through a packed bed of fixed spheres [27]. A turbulent regime with the standard k-model for the Re part numbers between 1000 and 2000 was considered in their study. The model was used to study fluid flow over a wide range of particle Reynolds numbers  for different packing arrangements. A unit cell (in red colour in Fig. 5) was composed of eight 1/8 spheres. Translational periodic boundary conditions were assumed with the pressure drop to be periodic at all unit cell faces. A semi-implicit method for a pressure-linked equation (SIMPLE) algorithm was used for correcting a pressure gradient based on the difference between the target mass flow rate and the actual one. Symmetry boundary conditions were imposed on all vertical faces bounding the fluid domain. The velocity and all gradients normal to the faces were set to zero. No-slip boundary conditions were applied on the surfaces of spheres [27].
In our simulation on lateral surfaces of the specimen, symmetry boundary conditions were chosen with the velocity normal to the surface, and all other normal gradients were set to zero. The velocity component parallel to the wall was computed only. In the inlet and outlet, periodic boundary conditions with a defined mass flow rate were chosen. On the sphere surfaces, no-slip conditions were defined where the fluid velocity was set to zero. The water entered the inlet at 70 ℃. The properties of water (with a proppant) were as follows: the dynamic viscosity 0.000406 Pa·s and the density 977.36 kg/m 3 . An unstructured tetrahedral mesh was also used to discretize the fluid domain. A very dense mesh was assumed in the vicinity of spheres, to capture high-velocity gradients and related high shear stresses. The mesh dependency test was carried out by measuring the pore average streamwise velocity Vy in the unit cell against the number of elements (Fig. 6). The streamwise velocity converged for the number of elements greater than one million (Fig. 6) and the mesh setting corresponding to one million elements was adopted for simulations.
Two numerical validation scenarios were considered on the cubic experiment specimen composed of eight 1/8-spheres with a diameter of 28 mm diameters (Fig. 5). The first scenario compared the calculated fluid velocities with the experimental [26] and the numerical results [27]. The second scenario compared the permeability results using the two most common equations: Carman-Kozeny (Eq. 12) and van der Hoef (Eq. 13). FVM was adopted to solve governing equations. The discrete system of linearized equations was solved using the algebraic multi-grid method by means of the incomplete lower upper (ILU) factorization technique. Hydrodynamic equations for u, v, w and p were solved as a single system by a coupled solver of ANSYS CFX [25]. A high-resolution scheme was used for the advection term. Pressure and velocity were coupled with a high-resolution scheme derived from the discretization algorithm developed by Rhie and Chow [47] and modified by Majumdar [48]. Figure 7a presents the dimensionless y-component of fluid velocities for Re part = 204.74 at y = 0 along the dimensionless x-value (blue line in Fig. 4). The distribution of the dimensionless y-component of velocity was parabolic and matched well with both the experimental [26] and numerical ones [27]. The dimensionless y-component of the fluid velocity at the center of the unit cell for five different values of the particle Reynolds number Re part is shown in Fig. 7b. The velocity variations are parabolic with increasing Re part . Good accordance of our numerical results with the experimental [26] and numerical ones [26] was obtained except for the result Re part = 12.17. Gunjal et al. [26] discussed in detail the possible causes of discrepancies between the results obtained in the experiment and those obtained in numerical simulations with Re part = 12.17. They came to the conclusion that the discrepancies were not related to numerical errors but to measurement inaccuracies in the experiment since a steady flow at a very low flow rate was difficult to maintain, and the greatest scattering of the experimental results was there for Re part = 12.17.
In the second validation scenario, the diameter of the beads at each corner of the unit cell ranged from 0.006 m to 0.038 m to get different porosities ϕ. The Darcy permeability coefficient K D was computed for each porosity from Eq. 6. To compute the Darcy permeability, the small pressure difference between inlet and outlet was imposed in such a way that the Re k < 0. (1 − ) 2 Figure 8 presents a comparison of the dimensionless permeability K D /D 2 according to our simulation (Eq. 6), the Carman-Kozeny equation (Eq. 12) and the van der Hoef equation (Eq. 13) for Re k < 0.1. Equation 13 is valid in the porosity range 0.35< < 0.99. To determine the average diameter for semi-spherical particles, the Sauter mean diameter [51] was assumed, defined as the diameter of a sphere of the same surface area/volume as a particle of interest. The calculated dimensionless permeability in Fig. 8 non-linearly increases with growing porosity. The numerical result is in good agreement with Eq. 12 in the range 0.30< < 0.77 and with Eq. 13 in the entire porosity validity range (Fig. 8). The differences between the simulation results (Eq. 6) and the results of Eq. 12 for ≤ 0.30 result from the overlapping of spheres whose shape was not perfectly spherical. For porosities higher than 0.85, Eq. 12 needs a correction factor.

Fluid flow in densely packed spheres
The validated full CFD model was compared with our 2D coupled DEM-CFD approach [21,22], based on a fluid flow network between granulates. The validation was performed under conditions typical for a hydraulic fracturing process. However, the hydraulic fractures were omitted in numerical simulations yet.

Specimen
The fluid flow characteristic was tested in a specimen consisting of 1.0 mm diameter spheres (Fig. 9). Some of the spheres overlapped one another. The designed specimen imitated a rock specimen. Discrete elements (spheres) were arranged in one layer in such a way that their centers of gravity were located in the specimen midplane. As a result, the specimen had one layer of dense packing of spheres with considerable overlaps of a maximum of 0.3 mm (Fig. 10b).
The size of the specimen was 10 × 10 × 1 mm 3 and the total fluid volume was 3.378e −8 m 3 . A small specimen was assumed to significantly shorten the computation time, particularly with the full 3D CFD model (see Sect. 5.4).
The specimen was composed of 131 spheres in a random arrangement. Other researchers who studied fluid flow in dense sphere systems, did not consider overlaps of spheres (e.g. [33,50,[54][55][56]) due to meshing problems. Very high pressures (up to 70 MPa) were chosen to reproduce the nature of a hydraulic fracturing process. The same sphere arrangement was used in both specimens. However, the specimen for simulations with the full CFD model solely contained the fluid domain with the fluid flowing around spheres that were initially used to create a model geometry of the fluid domain and then removed. The specimen for coupled DEM-CFD simulations was composed of both spheres (3D discrete domain) and fluid (2D continuous domain) (thus, two domains coexisted: the 3D discrete domain (sphere) and the 2D fluid domain located in the specimen mid-plane). The porosity was defined as the ratio of the fluid volume (void volume) to the total specimen volume (spheres volume and fluid volume) for both models. However, due to the two-dimensional nature of the fluid domain in the coupled DEM-CFD model, the fluid volume was differently computed. In the coupled DEM-CFD model, the spheres were interpreted as cylinders, i.e. as circles in the specimen mid-plane. The volume of cylinders is larger than that of spheres. Hence, the fluid volume in 2D was smaller than in 3D. This fact resulted in a different porosity in both models. The porosity of the DEM-CFD specimen was 13.864% while the porosity of the full CFD specimen was =34% and was higher than a typical one for rocks (e.g. the porosity of shales is typically less than 5%). Designing a lower porosity of the full CFD specimen was problematic regarding mesh generation (too many cells and too poor mesh quality). Even so, the fluid domain geometry was very complex due to overlapping spheres. Both the models could be compared for different porosities if the permeability was the same. Therefore, the coupled DEM-CFD model was calibrated to achieve the same permeability as the calculated one for the full CFD model. The 2D DEM-CFD model calibration procedure was based on the fact that in numerical simulations the height of the virtual channels (type S2S) was fitted to the target permeability. In spite of different Fig. 7 Dimensionless streamwise flow velocities V y /Vy in full CFD simulations as compared to experiment [26] (A) and numerical result in [24] (B): a for Re part = 204.74 at centre of unit cell (y = 0) along dimensionless x-value (R is the spheres radius) and b for 5 different Re part numbers at origin of unit cell

Initial and boundary conditions
The boundary conditions in the fluid domain of the full 3D CFD model are presented in Fig. 11. The turbulent, compressible fluid propagated through the rock specimen driven by a pressure drop between the inlet and outlet. The inlet and outlet boundary conditions were defined as 'opening, i.e. the fluid could flow through the boundary wall in both directions. The total inlet pressure was equal to 70 MPa and the average static pressure was defined at the outlet in such a way that the total outlet pressure was 30 MPa (the pressure difference was 40 MPa). In the outlet, the average static pressure over the entire outlet was specified in this way that the outlet pressure profile might vary but the average value was constrained to a specified value. At both the inlet and outlet, turbulence intensity defined as I = u /u (u -the rootmean-square of the turbulent velocity fluctuations and uthe mean velocity), was I = 5%. The turbulence intensity is a factor that shows velocity fluctuations and it can be used to determine a flow regime. Generally, a turbulence intensity of 1% or less corresponds to low intensity and greater than 10% corresponds to high intensity. At lateral surfaces of the specimen, free slip boundary conditions were defined with the velocity normal to the surface and wall shear stresses set to zero (i.e., v n,wall = wall = 0 ) to reproduce fluid flow close to unidirectional one at the macro-level and thus to facilitate the comparison. Only the velocity component parallel to the wall was computed. At the sphere surfaces and the domain corners, a no-slip condition was used. All the spheres were fixed and the water of 70℃ entered the inlet boundary. The material properties were as follows: the dynamic viscosity was μ = 0.000406 Pa·s and the reference density was 977.36 kg/m 3 . The reference pressure was set to 0.1 MPa in all simulations. The initial velocity was 0 m/s and the initial static pressure was 30 MPa in the fluid domain. The pressure difference of 40 MPa between the inlet and outlet resulted in a mass flow rate of 0.183 kg/s through the specimen. A fully implicit scheme was applied to a general discrete approximation of transient terms. Unsteady calculations were performed using a physical time step of Δt = 2·10 -6 s. The second-order backward Euler scheme was used for transient simulations. The convergence criterion at every time step was the maximum RMS value of the dimensionless velocity and pressure residuals or the maximum limit of 100 sub-iterations. The maximum RMS value was set to 1·10 -6 to ensure adequate damping of errors throughout the entire domain. The average values of about 100 subiterations were used to converge the solution per the physical time step. Each simulation was continued until a steady state was reached. In the DEM-CFD model, an adaptive method was employed to determine a physical time step that varied in the range from 2·10 -8 up to 2·10 -7 .
In the full CFD model, a very dense element mesh was used near surfaces of spheres to tackle high-velocity gradients in the boundary layer region. Usually, two main methods are applied for generating this type of meshes. For simple geometries like flow through spheres located at a certain distance from each other, thin layers of anisotropic prismatic or hexahedral elements near no-slip boundaries may be used to resolve sharp gradients in the wall-normal direction ( [57,58]). A coarser structured/ unstructured mesh may also be used in other domain regions far from surfaces of spheres to decrease the computational time. However, it is not an efficient meshing method in a dense and random packing of spheres located in an asymmetric arrangement with considerable overlaps. Those overlaps create several small regions and holes with acute angles between the adjacent spheres. For meshing this domain, a dense unstructured grid is used near all surfaces of spheres to efficiently cover boundary layers. Finally, sufficiently small unstructured linear tetrahedron elements were chosen for meshing the complex geometry. The element size was between 2.7·10 -5 m and 5.4·10 -5 m with a growth rate of 1.2. A close-up view of the mesh distribution on spheres' surfaces and the midplane is shown in Fig. 12. Small white spots on spheres in Fig. 12 are related to contacts between lateral slip walls and sphere surfaces (what may be seen in Fig. 11).
The same initial and boundary conditions for the fluid domain were defined in the coupled DEM-CFD model. The spheres in the DEM were fixed towards the axis OZ only. The cohesion assumed in DEM was high enough to significantly reduce the displacement of spheres. Therefore, the position of the spheres at the beginning and end of the simulations was almost the same.

Numerical results of full 3D CFD simulations
The commercial software ANSYS CFX [25] was solely used to simulate fluid flow in the continuous fluid domain. The spheres were initially applied for creating a model geometry of the fluid domain and next removed. For comparison purposes, the fluid flow contours in the full CFD model were presented in the specimen midplane. The change of calculation results in time was shown for 5 measurement points (Fig. 10b) whose coordinates were given in Table 2. A mesh dependency test for a turbulent flow regime of the full 3D CFD model was initially carried out by determining the pore average streamwise velocity Vy against the number of elements (Fig. 13). The velocity converged when the element number was greater than 13 million. An increase of the element number up to 15 million caused an effect of 1.0% only on the final average pore streamwise velocity. Thus, the mesh assembled of 15 million elements was used in this study. The transient simulation reached an equilibrium state after 0.07 ms in the full 3D CFD model, based on the mass flow rate at the outlet boundary (Fig. 14). Such a dense mesh was used to perform only one simulation. It is obvious that simulation of a varying geometry process like hydraulic fracturing with commonly used numerical methods such as FVM, LBM, SPH, etc. are not almost feasible due to the extremely long computational time related to a remeshing process. Hence, the model of the fluid flow network (e.g. [11,13,21,22]) constitutes an effective method to overcome those time problems.
The 3D streamlines representing the velocity in the full CFD model are shown in Fig. 15. The flow accelerated at the upstream surfaces of the spheres due to the reduction of the cross-sectional area in the main flow direction and slowed down at the downstream surfaces of the spheres. The interaction of streamlines in some regions caused a generation of swirling flow (e.g., bottom right corner of the specimen Fig. 15). In those regions, the flow direction changed and rotated around the axis along which it moved under the influence of angular momentum. The streamlines were directed in a negative y-direction in some regions. Figure 16 shows the total pressure distribution in the 3D domain using the full CFD model. The boundary condition imposed the pressure of 70 MPa in the inlet and 30 MPa in the outlet. The average total pressure over the entire fluid domain was 50 MPa. The pressure increased at the front and both sides of the sphere bounds and decreased at the backside. The presentation of contours on a 2D mid-plane was obscure due to the high porosity (Fig. 10b). Therefore, to compare more precisely the results, the 2D contours were plotted in Statistica with the spheres were removed from the plot. In addition, the variability of the variables was smoothed using the distance-weighted LS method to obtain consistent isolines.

Comparison of numerical results between full CFD model and DEM-CFD model
In general, in the DEM-CFD model, bonds between particles may be broken during simulations (Appendix 'A'). They were not broken in a deforming non-homogeneous granular specimen considered here since fractures were not taken into account yet (Sect. 1). However, the contact stresses changed during simulations that resulted in a change of the channel height in the fluid flow network and in a change of virtual pore volumes. Consequently, the fluid flow characteristics varied during the simulation until it reached an equilibrium state (Sect. 3). Figure 17 compares the pressure distribution at the midplane in the full CFD model (Fig. 17a) with the pressure contours obtained in simulations with the coupled DEM-CFD model (Fig. 17b) after reaching the equilibrium state. Both models provide the same pressure distribution. Due to a non-symmetric arrangement of spheres, the pressure isolines, especially in Fig. 17b, are not straight horizontal lines. The isolines in Fig. 17a are similar in the entire domain while in Fig. 17b they are close to each other near Fig. 12 View on unstructured tetrahedral mesh (with zoom) assumed in granular rock specimen the inlet (and outlet) and the distance between them gradually increases, particularly at the left side of the specimen. The overall flow resistance due to shear stresses remains at the same level in the full CFD model while it is higher in the region near the inlet and outlet in the DEM-CFD model. It should be noted that the fluid flow network in the DEM-CFD model is not orthogonal. Significant discrepancies in the pressure distribution may, thus, indicate shortcomings in DEM-CFD models.
The evolution of the total pressure in time is presented in Fig. 18 for 5 points of Fig. 10b. There exists a good agreement between the final steady-state results at those 5 points. The full CFD model simulation reached the equilibrium state after 0.07 ms. However, the coupled DEM-CFD simulation reached this state after 0.00025 ms. The water volume in the fluid domain for the full CFD model was 3.478e −8 m 3 while the water volume for the fluid flow network model was solely 6.89e −12 m 3 . This difference in the water volume strongly influenced the inertia forces in the DEM-CFD model and consequently the fluid velocity. The fluid flow in channels in the DEM-CFD model was 2D laminar without disturbances whereas in the full CFD model, it was the 3D flow with vortices and disturbances. Therefore, the CFD-DEM model reached the equilibrium state 280 times faster. Figure 19 demonstrates the density distribution in the full CFD model that is similar to the pressure distribution (Fig. 16) since the density solely varies with pressure variations in our barotropic model. Near the inlet where the maximum pressure and consequently maximum density change existed, the density solely grew by 2.85% as compared to the reference density of 977.36 kg/m 3 .
To investigate the accuracy of the barotropic equation of state, another simulation of the 131 spheres sample with the full CFD model was performed based on the IAPWS material library. Table 1 shows the density variation, pore average density, streamwise velocity on the outlet boundary, and pore average streamwise velocity for both the barotropic and IAPWS model. In addition, the results for an incompressible model were also attached. There are small differences between the range of density variations and pore average densities in the two compressible models. Both the models share the same final Vy at the outlet as well as Vy . The values are almost the same as compared with an incompressible model. Hence, the density variations are not significant and the flow can be described using the incompressible fluid model in a considered pressure range.
The density distribution at the mid-plane of the full CFD model (Fig. 20a) was compared with the density distribution of the coupled DEM-CFD approach (Fig. 20b). Both the models resulted in the same density distribution which corresponded to the 2D pressure distribution (Fig. 17). The only difference is in the density isoline location (Fig. 20). They are far away from each other as compared with the pressure isolines (Fig. 17). The discrepancy in the pattern of the fluid density isolines results from a greater concentration of discrete elements on the right side of the specimen. Through calibrating the DEM-CFD model, the mass flow rate in this region is almost the same for both models. Consequently, the pressure distribution is also almost the same. However, the fluid velocity is significantly different that can lead to small density differences between the models. The density variation in time in 5 measurement points (Fig. 10b) is demonstrated in Fig. 21 for two approaches: full CFD and simplified DEM-CFD. The density variation is similar at all points independently of the approach.
The 3D streamwise velocity distribution based on the calculation results of the full CFD model is demonstrated in Fig. 22. High fluid velocities occurred close to lateral walls where the slip boundary condition was applied, reaching a maximum value of 205 m/s. The negative velocities were observed in some areas in the contour plot. The pressure difference between the inlet and outlet resulted in the average streamwise velocity of 20.65 m/s at the inlet/outlet boundary surface. Usually, the velocity inside porous media is presented in the form of superficial velocity U s . This velocity which is the average value of the streamwise velocity in the total specimen volume (i.e., fluid and spheres) showed a constant velocity all around the domain (U s was 18.12 m/s in the granular rock specimen).
The Re k and Re part numbers in the full 3D CFD model were calculated as 1480 and 15,784 for the average volume density of 996 kg/m 3 . Due to overlaps, the final shape of the particles was not spherical and the Sauter mean diameter [53] was used to obtain an average diameter for all spheres. The Sauter mean diameter of all 131 spheres was 1.22 mm and this value was used to evaluate the value of Re part . Since the fluid in the coupled DEM-CFD model flew in narrow channels, the velocity range was totally different as compared with the full CFD model. The fluid moved much faster along the channels of the fluid flow network to achieve the same mass flow rate as in the full CFD model. The In the 2D DEM-CFD model, the maximum fluid velocity (1489 m/s) was 10 times higher than in the 3D full CFD model (153 m/s) and the Reynolds number was 10 times lower (1671) than in the full CFD model (15,784). The fluid flow still remained in a laminar flow regime while the fluid flow was turbulent in the full CFD model. This difference is due to a different formulation of two models; the full CFD model realistically reproduced the fluid domain geometry while the fluid domain was simplified in the DEM-CFD model by applying a fluid flow network composed of simple low-height channels. This resulted in a difference in the fluid volume and specimen porosity. The specimen porosity was 34% (full CFD) and 13.864% (DEM-CFD). The specimen porosity in the 3D full CFD model was almost the smallest one due to mesh generation problems while the specimen porosity in the DEM-CFD model was the highest to preserve the laminar flow. Increasing the specimen porosity in the DEM-CFD model may lead to the Reynolds number exceeding the laminar flow limit and thus contributing to the model incorrectness (the model calibration is limited by the critical Reynolds number for laminar flow in channels). In the DEM-CFD model, the fluid velocity, Reynolds number and permeability are strongly affected by the channels' height. Therefore, the channels' height is used to calibrate the DEM-CFD model for the assumed permeability. Consequently, increasing the assumed permeability (for the greater specimen porosity) requires a growth of the channels' height that causes an increase in the fluid velocity. To overcome this shortcoming of the DEM-CFD model, the fluid flow model should be modified to be closer to the full CFD model.
The Reynolds number Re in the DEM-CFD approach was computed in the channels using equation Re ch = V ch h , Even though the fluid velocity was much higher in the DEM-CFD approach, the Reynolds number was lower since the channel height was very small (4.57·10 -7 m). The calculated pressures and densities at 5 measuring points in Fig. 10b are summarized in Table 2 for the full CFD approach and the coupled DEM-CFD approach.
The vorticity is the quantity that shows the rotation of a fluid flow direction. Figure 23 presents the 3D out-of-plane streamwise vorticity distribution in the full CFD model with a clockwise rotation corresponding to positive vorticity. Some eddies in a 2D vertical mid-plane are shown in Fig. 24. Alternating right and left rotation regions appeared at the left and right sides of each sphere. The regions of concentrated vortices were located in the shear (boundary) layer of each sphere. Interaction between the shear layers generated eddies of various sizes in the specimen.
For very low pressures, the flow disturbances and turbulent kinetic energy immediately dissipated and flow remained laminar. The laminar model was used for very low pressures to find Darcy's permeability in a creeping flow regime. To evaluate permeability using Darcy's law, the pressure drop between the inlet and outlet was reduced   The fluid flow characteristics depend on two important fluid properties: viscosity and density. In turbulent flow, it also depends on other parameters, depending on the turbulence model used. The other calculations showed that a viscosity increase caused a slight reduction in the pore average streamwise velocity. The density change had a more pronounced effect on the pore average streamwise velocity than the dynamic viscosity. In addition, the slope of density changes was sharper for smaller densities.
Comparing the performance of both approaches, for the full CFD model, the time to generate the mesh was about 1 h using 18 computer cores and the computation time was about 60 h using 144 computer cores. In the simplified coupled DEM-CFD approach, the total time for the automatic fluid flow network generation and calculations was about 52 h using 1 computer core. Thus, the coupled simplified DEM-CFD simulations are feasible with regard to time in contrast to full CFD calculations.

Conclusions
Single-phase fluid flow through granular body composed of fixed densely packed overlapping spheres imitating a non-homogeneous rock specimen under high pressures in isothermal conditions was numerically analyzed using the 3D full CFD model in a continuous domain between spheres. The numerical CFD results were compared with the corresponding results of a simplified coupled 2D DEM-CFD approach using a virtual pore network in a domain between spheres that was composed of channels. The following conclusions may be offered: • Both the numerical approaches provided almost the same pressures, densities and mass flow rates. However, the maximum flow velocity in the DEM-CFD model was about 10 times higher than in the full 3D CFD model to transport enough mass for achieving the same mass flow rate. Velocity differences were due to different water amounts in the 3D specimen (water flow in the entire continuous domain between spheres) and 2D (water flow in the fluid flow network between spheres only). Thus, the simplified coupled DEM-CFD approach may be used in hydro-mechanical simulations. • At the considered pressure difference (40 MPa) and temperature (70℃), the changes in density were insignificant and water might be assumed to be incompressible in both models. • Due to the lack of experimental tests for high-pressure flow at the meso-scale, the presented procedure in Sect. 5 may be applied to validate indirectly any 2D/3D DEM-CFD models, based on the concept of a fluid flow network. • The performance of the 2D DEM-CFD model was significantly faster (52 h using 1 core) as compared to the full 3D CFD model in a continuous domain between spheres (60 h on 144 cores).

Appendix 'A'
In DEM, particles interact with each other during translational and rotational motions through a contact law and Newton's 2nd law of motion using an explicit time-stepping scheme. A cohesive bond is assumed at the grain contact The DEM model for cohesive-frictional materials is described in detail in [21,22,34,35]. It can be summarized as follows (Eqs. A1-A7, Fig. 25) : F n × tan c ≤ 0 (after contact breakage) where �⃗ F n -the normal contact force, U-the overlap between discrete elements, �� ⃗ N-the unit normal vector at each contact point, �⃗ F s -the tangential contact force, �⃗ F s,prev -the tangential contact force from the previous iteration, �⃗ X s -the relative tangential displacement increment, K n -the normal contact stiffness, K s -the tangential contact stiffness, E c -Young's modulus of the particle contact, ν c -the Poisson's ratio of particle contact, R-the particle radius, R A and R B -the contacting particle radii μ c -the Coulomb interparticle friction angle, F s max -the critical cohesive contact force, F n min -the minimum tensile force, C c -the cohesive contact stress (maximum shear stress at pressure equal to zero) and T c -the tensile normal contact stress, �⃗ F k damp -the damped contact force, �⃗ F k and � ⃗ v k p -the kth components of the residual force and translational particle velocity v p and α d -the positive damping coefficient smaller than 1 (sgn(•) returns the sign of the kth component of velocity). Five main local material parameters are required for DEM simulations: E c , υ c , μ c , C c and T c . In addition, the particle radii R, particle mass density ρ c and damping parameters α d have to be known. In general, the material constants are identified in DEM with the aid of simple laboratory tests on the material (uniaxial compression, uniaxial tension, shear, biaxial compression) [21,22,34]. The calculations were carried out (A6) F s max = C c R 2 and F n min = T c R 2   [23,24] the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.