Coupled flow and deformation processes are important for a variety of applications spanning the multiphase multicomponent flow models in reservoir simulation [1][2][3][4][5], through subsidence and consolidation in soils [6,7], to the phenomena in human physiology including eye modeling [8], modeling bone fractures and prosthetics [9], and poro-elastic brain models [10]. We refer, e.g., to the recent collection [11] featuring extensive references to the rich variety of applications and methods.
These phenomena received considerable attention from computational scientists and modellers, and there is substantial work devoted to the construction and analysis of numerical schemes as well as their analyses; see, e.g., the work in [4,10,12,13]. While the construction or robust schemes is important, the predictive power of particular simulations depends primarily on the ability of a model to represent a physical phenomenon and the availability of realistic data. Furthermore, the typical time and spatial scales of such processes vary widely, which should be reflected in the model assumptions and data. In general, these time and spatial scales are harder to predict for coupled systems than for individual process, and in particular are harder to find for the coupled flow and deformation than for the flow alone. In this paper we propose a clear novel workflow to address the connection between the pore-and Darcy-scales for flow and deformation. Our algorithms derive from the combined expertise in granular mechanics, computational fluid dynamics (CFD), and upscaling.
Darcy scale flow and deformation model. In particular, in this paper we consider an extension of the well known linear Biot model solved for the displacement 𝑢, pressure 𝑝, and flux 𝑞 𝑓 in a Darcy scale domain 𝛺 ⊂ R 𝑑 -∇ ⋅ [𝜆(∇ ⋅ 𝑢)𝐼 + 2𝜇𝜖(𝑢)] + 𝛼∇𝑝 = 𝑓 + ρ𝐺∇𝐷, (1a) Sand (medium/ fine): 0.25 [14](Pg. 407) Silt : 0.30-0.35 [14](Pg. 407) Clay(medium/firm): 0.30 [14] Sand: 10 -13 -10 -11 [m 2 ] (saturated) [14](Pg. 373) Silt: 10 -15 -10 -13 [m 2 ] (saturated) [14](Pg. 373) Clay: 10 -18 -10 -15 [m 2 ] (saturated) [14] 𝜕 𝑡 (𝑐 0 𝑝 + 𝛼∇ ⋅ 𝑢) + ∇ ⋅ 𝑞 𝑓 = 𝑓 𝑓 , (1b)
complemented by appropriate initial and boundary conditions. The data for the problem include the Lamé coefficients 𝜆, 𝜇, the specific storage coefficient 𝑐 0 which involves the porosity 𝜂, the fluid viscosity 𝜇 𝑓 and density 𝜌 𝑓 , and permeability 𝜅. See Table 1 for notation. Typically the data 𝑐 0 , 𝜂, 𝜅 for modeling efforts of (1) at the Darcy scale are calibrated from experiments. Data at Darcy scale from upscaling. Much work has been devoted to getting the Darcy scale model data 𝜂, 𝜅 via upscaling from the pore-scale; see, e.g., well known homogenization efforts for pore-to-Darcy in [18][19][20][21] as well as our own work involving xray micro-CT images in [22][23][24][25][26]. The imaging of some porous medium volume 𝜔 provides a grid with voxel resolution 𝑀 𝑥 × 𝑀 𝑦 × 𝑀 𝑧 of the image of 𝜔, with each voxel assigned 1 or 0 depending on whether it is recognized as part of the rock 𝜔 (𝑟) and fluid filled 𝜔 (𝑓 ) subdomains, with the porosity 𝜂 = |𝜔 (𝑓 ) | |𝜔| . The workflow for such efforts typically involves the following steps Imaging 𝜔 (2a)
→ pre-processing (2b)
→ upscaling → (𝜂, 𝜅). (2d) Such efforts are especially useful for the understanding of the changes in permeability when the pore-scale geometry is changing due to some coupled process such as precipitation, microbial plugging, or phase transitions. The computational effort in ( 2) is typically dominated by the time spent on imaging the Representative Elementary Volume (REV) 𝜔 in (2a) and by the computational complexity of CFD in (2c), which typically scales as 𝑂(𝑀 2 ) due to the cost of a linear solver. (Here 𝑀 = 𝜂𝑀 𝑥 𝑀 𝑦 𝑀 𝑧 is the number of fluid cells of the computational grid covering the flow portion 𝜔 (𝑓 ) of 𝜔.) This high computational complexity suggests to keep the resolution of the grid covering the REV and/or the size of REV small. At the same time, the accuracy of the process of upscaling [17,19] calls for the REV to be large enough to faithfully represent the average qualitative and quantitative processes, while the accuracy of the CFD simulations calls for large enough 𝑀. In the end, we work in practice with moderate 𝑀 = 𝑂 (10 6 ).
In this paper we are interested in the modifications of the pore-scale domain 𝜔 due to the deformation and the associated change of the permeability 𝜅 𝜅 = 𝜅 ( deformation ),
while keeping other data for (1) including 𝜆, 𝜇, 𝑐 0 fixed. In fact, in a saturated porous geologic material the bulk stiffnesses of the fluid and the grains are many orders of magnitude greater than the bulk stiffness of the granular skeleton. Thus, changes in permeability are dominated by rearrangement of the pore network in response to the deformation of the solid skeleton. A common strategy to obtain (3) is to predict local changes in porosity 𝜂 due to deformation, and to infer (3) from the well-known Carman-Kozeny relationship
in which 𝛼 𝐶𝐾 is a positive parameter that needs to be determined, and 𝜅 0 is the reference conductivity. In practice, the actual values of 𝜅 0 , 𝛼 𝐶𝐾 are assumed. Once ( 3) is known, one can use it in a nonlinear extension of (1). There has been work on this particular direction, primarily in the context of visco-elasticty in human tissue [8,27]. The specific form of (3) is typically postulated, e.g., with (4), rather than derived. The exception is the work in [28][29][30] where the relationships are defined from the pore-scale problems under periodicity assumptions. While these are useful from a theoretical standpoint, it is not easy to get the data for (3) from these pore-to Darcy-scale efforts for the specific scenarios of interest.
In turn, in principle, a relationship such as (4) and in particular the value of 𝛼 𝐶𝐾 could be found following the process (2) based on a sufficiently rich collection of images of some REV 𝜔 under deformation. However, such imaging efforts under deformation are usually not feasible. In this paper we propose instead to derive (3) or ( 4) for the Biot model (1) using the discrete element method (DEM). The algorithms connecting DEM to (3) are the focus of this paper and the presentation of the practical algorithms and their tests and results is our major goal.
Discrete element method (DEM) has been used successfully for a variety of predictive models to produce realistic assemblies of porous grains, i.e, different REV 𝜔, in response to an external load, under various assumptions. DEM models are used to simulate the emergent mechanical behavior of granular materials (e.g., sands, grains, powders) and bonded particulates (e.g., concrete, rock). In a DEM simulation, the response of a particulate assembly to a physical excitation (body forces and boundary tractions) is predicted based upon simultaneous solution of Newton's equations of motion for every particle in the assembly subject to well-defined interaction laws at particle contacts. That is, given an ensemble of 𝑁 𝑝 interacting particles, a general computational approach consists of a repetition of the following steps: (a) identification of the locations of all particle-particle contacts and computation of the corresponding normal vector orientations; (b) application of the contact force-displacement law for determination of the interaction forces between particles; and (c) time integration of the equations of motion to determine updated particle state. DEM models are particularly well-suited for simulation of emergent phenomena that elude ready analytical or constitutive description, such as the equilibration of a collection of spheres under gravity within a container. DEM is based on compatibility of particulate displacements through simple contact relationships rather than an overarching constitutive law -only the basic interaction physics and the compatibility criteria must be satisfied for a given time step.
Objectives of this work. In this paper we propose a ''proof-of-concept'' approach connecting the DEM to the upscaling efforts in order to get (3). We do not aim to work with any particular granular material or soil. However, we employ realistic data and realistic deformation scenarios.
In our approach we follow (2) but modify step (2a) to involve DEM instead of imaging:
Specifically, in (5a), we consider spherical particles of various sizes which are not individually deformable; we assume that the particles are governed by Newton's laws, and we use DEM to simulate the deformation of the solid phase corresponding to some assumed deformation scenario. The set of particles forms the solid phase 𝜔 (𝑟) which is found in step (5b), from which we also get 𝜂. Next, in step (5c) we use the well-known Marker and Cell (MAC) scheme to simulate the Stokes flow in 𝜔 (𝑓 ) . We obtain the conductivity 𝐾 = 𝜇 -1 𝑓 𝜅 by upscaling in step (5d) the fluid velocity and pressure found from step (5c). Actually, to make our results meaningful, we repeat the entire process several times in response to some scenario of deformation, and in the end we obtain a collection of DEM-derived assemblies (𝜔 𝑘 ) 𝑘 , with each 𝜔 𝑘 composed of some 𝑁 𝑝 particles responding to some given load, e.g., increasing the load between step 𝑘 and 𝑘 + 1. Typically, with increased load, the change from 𝜔 𝑘 → 𝜔 𝑘+1 is due to either individual particles of the rock phase ''settling'', the rock phase undergoing deformation, or some combination of both.
In the end for each 𝜔 𝑘 we obtain the corresponding (𝜂 𝑘 , 𝜅 𝑘 ). We can use their collection to calibrate some parametric model for (3) e.g., (4), and carry out the simulations of Biot model (1) which we extend to allow the dependence of permeability on deformation.
In summary, the workflow (5) we propose has various advantages over other approaches including experimental. Collection of real 3D pore geometry information is expensive and, for specimens constructed in the laboratory, still not guaranteed to reproduce the intricacies of the pore structure that exists in-situ. Conversely, it is not always possible to demonstrate the physical or thermodynamic admissibility of synthetic pore structures generated via artificial means. DEM simulations can bridge this gap by computationally generating pore structures that are both physically realistic and have an appropriate level of geometric and topological complexity. While the computational costs are non-negligible (hours to days for a single assembly of tens of thousands of particles for DEM alone) followed by hours to days of CFD and upscaling, they can be performed largely unattended on very inexpensive commodity workstations.
Outline. In Section 2 we describe DEM and its use for (5a). In Section 3 we describe the pre-processing step (5b). Section 4 is devoted to (5c) and solving the flow model, with particular attention paid to the workflow of obtaining the solutions on the geometries 𝜔 generated by DEM. In Section 5 we outline the method of upscaling in step (5d) and provide the results. In Section 6 we use this data for Biot simulations. We conclude and summarize in Section 7. Appendix contains supporting data.
Remark 1.1 (Notation). We adopt the following notation. For a set 𝐴 we denote by |𝐴| its measure, and by Ā its closure: we annotate the difference between 𝐴 and Ā but occasionally when taking union or set differences we abuse the notation rather than annotate the specific differences. In particular, the domains for PDEs such as 𝛺 in Biot model (1) or the flow domain 𝜔 (𝑓 ) for the Stokes flow model are open. However, to avoid clutter, we do not distinguish these from the computational domains, a union of computational cells which are closed.
We denote by 𝑁 𝑝 the number of particles 𝐵 𝑝 , 𝑝 = 1, 2, … , 𝑁 𝑝 modeled by DEM which are closed balls whose union forms 𝜔 (𝑟) .
In this Section we describe step (5a) of creating realistic synthetic porous media microstructures using DEM simulations. We start with literature overview, and next proceed to define the steps towards (5a).
DEM has been used to simulate particulate assemblies in two (e.g., [31,32]) and three (e.g., [33][34][35][36]) dimensions, using spheres, disks, ellipses, ellipsoids, and clumps to represent individual particles. A large number of physical processes have been simulated using the DEM, including particle breakage (e.g., [37,38]), rock fracture [39,40], anisotropic continua [41], and fruit damage due to impact [42]. One area in which the DEM has been heavily exploited is in the study of discrete material microstructure [33,43]. Because many parameters may be tracked with DEM that are not readily observable in the laboratory (e.g., particle orientation, force chains), the DEM is an excellent tool for studying and quantifying the microstructure of discontinuous systems. Previous studies have quantified the effects of particle shape on fabric evolution [44], the use of the anisotropy of contact orientations to predict large-strain shear strength [45], and particle rotations [46]. The DEM has been experimentally validated for the two-dimensional problem of fabric evolution in long rods subjected to boundary forces [47] and in three dimensions for regularly packed assemblies of steel balls [48]. The DEM has also been used in conjunction with morphological and stereological methods to successfully quantify evolving microstructure in two dimensional particulate assemblies [49,50].
The global response of a particulate assembly in DEM is allowed to develop naturally based on element interactions, resulting in the organic simulation of emergent phenomena. This is a significant advantage over continuum models (e.g., Finite Element Method) for simulating problems with very large strains or discontinuities (e.g., fractures) because there is no numerical instability or need for re-meshing in these cases. The number 𝑁 𝑝 of particles considered in a simulation generally ranges from several thousand to several tens of thousands but might be as high as tens of millions. There is no actual upper bound on the number of particles as memory demands are modest; the limit is practical because computations are processor-bound and the cost scales steeper than linearly with the number of particles.
Model inputs for DEM are largely comprised of well-known or readily accessible material properties for individual particles within the assembly. For example, a sand grain might have the properties of pure quartz: specific gravity of 𝐺 𝑠 = 2.65, shear modulus of 𝐺 = 29 [GPa], Poisson ratio of 𝜈 = 0.31, and a coefficient of sliding friction of 𝜙 = 0.5. Remaining inputs include model geometry, particle size distribution, and boundary conditions, which can be defined to describe the problem of interest.
Assemblies of particles are generally equilibrated within the model domain subjected to a body force and then exposed to an external stimulus (e.g., boundary tractions). During the simulation, the location of every particle and every contact is known, as are the force, normal vector, and shear vector at every contact. This information can be combined and post-processed to quantify assembly fabric (microstructure) and stress state (local and global) -and, importantly, their evolution with time.
Z. Hilliard et al.
We describe now our use of DEM, assuming some specific parameters. The choices of parameters will be revisited in Section 2.3. We select a model size (i.e., volume) and particle gradation that will produce a simulation having a number of particles 𝑁 𝑝 ≈ 12, 000 that is reasonably tractable for steps (5b) and (5c) of our workflow.
In this paper we consider spherical particles 𝐵 𝑝 , 𝑝 = 1, … , 𝑁 𝑝 : each particle is defined by its radius 𝑟 𝑝 and position 𝑥 𝑝 . In each DEM simulation, a particle assembly is created by placing the individual spheres into a bounding cube X until a desired porosity is reached. We draw particle radii from a log-linear distribution given by:
where 𝑟 𝑚𝑖𝑛 and 𝑟 𝑚𝑎𝑥 are minimum and maximum radii and is the uniform distribution.
Similarly, the three Cartesian components of the particle centroids (𝑥 𝑝 ) 𝑝 are independently drawn from a uniform distribution over X with dimensions reduced by the mean particle radius at each face of the box to prevent particle escape. In this approach, particles are ''hard core'' random: no measures are taken to prevent particle overlap, but no two centroids may occupy the same point in space. Since repulsive energy between spheres is proportional to overlap, significant contact forces are generated. The system is gently relaxed by stepping forward in time while zeroing the kinetic energy of all particles every few time steps. Once sufficient energy has been dissipated, the system is allowed to progress to equilibrium normally.
Once in equilibrium, the assembly has (nearly) exactly the specified porosity, but some non-prespecified stress state that is a function of the size distribution, particle mechanical properties, and contact force network topology. A numerical servo algorithm is then used to achieve a hydrostatic state of stress (𝜎 0 = 50 [kPa] in our case). This numerical servomechanism constantly adjusts wall velocity as a function of the difference between the actual and desired boundary stresses, while ensuring that the velocity remains below the value required for stability in any given timestep.
We note here that it is effectively not possible to create an assembly of particles that simultaneously possesses a specified gradation, specified porosity, and specified stress state -at least one of these requirements must be relaxed. Nonetheless, the change of porosity associated with achieving the desired stress state is small, typically < 1%. The system is now prepared for the ensuing perturbations. We denote the location of the centroid of each particle at this state 𝑥 0 𝑝 . Starting from this initial equilibrium, we compress the system vertically under zero lateral strain conditions. The same servo algorithm mentioned previously is used to gradually achieve a vertical stress of 𝜎 𝑣 = 100 [kPa] (±0.5%). The system is again cycled to kinetic equilibrium, and the vertical stress is measured. If the equilibrated stress is more than ±0.5% from the specified value, the position of the top face of the wall is adjusted until once again reaching the desired state of stress. This cycle is repeated until the desired stress state is maintained once the system has reached kinetic equilibrium. The coordinates of all particle centroids (𝑥 𝑘 𝑝 ) 𝑝 are recorded at this point and the system is then taken to the next point of desired equilibrium. Specifically, we capture snapshots of particle locations every 100 [kPa] up to 1000 [kPa] and then similarly back down to 100 [kPa] denoting each of these steps by numbers 𝑘 = 1, 2, …. All simulations are performed in the absence of gravity. In other words, we have generated the snapshots ((𝑥 𝑘 𝑝 )
We use this same approach to explore six scenarios of assembly response of an assembly of spheres to applied boundary tractions under different assumptions. We describe these scenarios below numbered (A), (B1), (B2), (C), (D), (D.r).
In this simulation, we start with an initial configuration of spheres in the bounding cube and simulate the bounding box as being comprised of six perfectly smooth (i.e., frictionless) walls. We consider this to be the ''baseline'' simulation.
This simulation is the same as Case (A) except that a greater number of calculation cycles are allowed as the system seeks kinetic equilibrium at each load level. This scenario is intended to confirm that the baseline state of equilibrium is sufficient for the ensuing calculations.
This simulation is the same as Case (B1) except that a different value is used to seed the random number generator at the beginning of the calculation. In this way, we confirm that system response is not unduly influenced by the values of the random variables so long as the same generation and equilibration procedures are used for every simulation.
This simulation is the same as Case (A) except that the particle volumes are cut in half (i.e., 𝑟 𝐶 = 𝑟 𝐴 ∕2 1∕3 ), so this case has 𝑁 𝑝 ≈ 24, 000 spheres instead of 𝑁 𝑝 ≈ 12, 000 in Case (A). We use this scenario to confirm that the results of the simulations do not depend on the size of the model, and that Case (A) has a significant enough number of spheres to constitute an REV.
This simulation is the same as Case (A), except the walls used to form the bounding box are frictional. Frictional walls serve to ''lock in'' stresses during loading such that maximal hysteresis is observed during unloading.
Z. Hilliard et al.
This case is not a new simulation but a further exploration of the results from Case (D). We take the snapshots from (D) and remove spheres to increase the porosity of the synthetic porous media. For the snapshots at 100 [kPa], 500 [kPa], and 1000 [kPa] (all taken from the sequence where the load is increasing), four sets of spheres were removed.
With the spheres in (D) labeled 𝐵 𝑝 , 𝑝 = 1, 2, 3, … , 12000, the first set of spheres to be removed corresponds to indices 𝛱 (1) 𝑟𝑒𝑚𝑜𝑣𝑒𝑑 = {𝑝 (1) 1 , 𝑝 (1) 2 , 𝑝 (1) 3 , … , 𝑝 (1) 𝑁 (1) 1
}, the second to 𝛱 (2) 𝑟𝑒𝑚𝑜𝑣𝑒𝑑 = {𝑝 (2) 1 , 𝑝 (2) 2 , 𝑝 (2) 3 , … , 𝑝 (2)
(2) 2
} and so on where 𝑁 (1) 1 < 𝑁 (2) 2 < 𝑁 (3) 3 < 𝑁 (4) 4 . These are chosen so that the set of spheres removed is consistent across snapshots with 𝛱 (1) 𝑟𝑒𝑚𝑜𝑣𝑒𝑑 ⊂ 𝛱 (2) 𝑟𝑒𝑚𝑜𝑣𝑒𝑑 ⊂ 𝛱 (3) 𝑟𝑒𝑚𝑜𝑣𝑒𝑑 ⊂ 𝛱 (4) 𝑟𝑒𝑚𝑜𝑣𝑒𝑑 . In the end, the scenarios listed above give each the data ((𝑥 𝑘 𝑝 )
) 𝑘 which constitutes complete geometrical information about the corresponding REV 𝜔 𝑘 .
We comment now on the specific choices of DEM parameters we made. These are based on our experience with a variety of scenarios from which we selected those discussed above (A), (B1), (B2), (C), (D), (D.r).
First, particle-level mechanical properties such as 𝜙 influence the packing of the particles in the assembly and, consequently, the contact force network in the solid phase of the system. Selection of different parameters would lead to different global porosities with different local geometries, but they would be no more or less ''correct'' than those we used. That is, both are physically admissible, they are just representative of different granular materials. In the end, then, we elected to use properties consistent with a silica sand.
Next, the number of particles was selected to provide a balance between the required complexity and computational tractability. The reasonableness of our selection may be justified through consideration of Case (A) with 𝑁 𝑝 = 12000, and Case (C) with 𝑁 𝑝 = 24000. Flow results from these two cases are nearly the same, implying that 𝑁 𝑝 = 12000 is an appropriate REV.
The initial mean stress (i.e., 𝜎 0 = 50[kPa]) was selected to be smaller than the initial stress 𝜎 1 = 100[kPa] used for flow calculations to avoid any stress history effects in the assemblies. Specifically, since all assemblies approach their first measurement point of 100 [kPa] from an initial equilibrium of 50 [kPa], we can ensure that any stress effects associated with initial generation are appropriately erased (see also Case (B2) relative to Case (B1).
The number 𝐾 = 19 and the range of loads for the snapshots was selected to be generally representative of near-surface geologic stress conditions, corresponding to geostatic stresses at depths from approximately 10 to 100[m].
In this Section we describe how to conduct the step (5b) following (5a). After (5a) we know the position of all particles {{𝐵 𝑝 } 𝑁 𝑝 𝑝=1 } and we can prepare the computational grid for the flow solution in 𝜔 (𝑓 ) . By design 𝜔 = 𝜔 (𝑓 ) ∪ 𝜔 (𝑟) , and clearly 𝜔 (𝑟) should represent ⋃ 𝑝 𝐵 𝑝 . However, there are two practical concerns. First, DEM uses a soft-contact approach to estimate the inter-particle forces. This introduces a slight overlap between some spheres so that |𝜔 (𝑟)
| < ∑ 𝑝 |𝐵 𝑝 |. This overlap is not significant in practice: in our simulations, we have consistently approximately only 0.01% of the domain occupied by overlapping spheres.
Second, for flow simulations we choose hexahedral grid over
where each 𝜔 𝑖𝑗𝑘 is a closed cube. This choice is natural for the workflow (2) and for the pore-scale domains obtained from imaging where the data is provided in the form of voxels, i.e., binary values indicating for each voxel ω𝑖𝑗𝑘 whether it belongs to 𝜔 (𝑓 ) or 𝜔 (𝑟) .
In this paper we follow this choice also for ( 5) even though we cannot represent any sphere exactly on a hexahedral grid, and have to work with approximations to 𝐵 𝑝 . This choice is justified based on our prior experience for (2) working with synthetic data or with spherical glass beads data forming 𝜔 (𝑟) , since with adequate resolution we encountered only negligible difference between the use of hexahedral vs body-fitted grids, while the 𝑖𝑗𝑘 numbering and the simplicity of flow algorithms had significant advantages of the former. Thus in the end we work with hexahedral grid approximations 𝜔 (𝑟,′) to 𝜔 (𝑟) . In fact, an additional step of approximations 𝜔 (𝑟,′) → 𝜔 (𝑟,′′) is needed to make the flow simulations robust in the corresponding 𝜔 (𝑓 ,′′) ; clearly, we attempt to make 𝜔 (𝑟,′′) as close as possible to the ''true'' 𝜔 (𝑟) defined by ⋃ 𝑝 𝐵 𝑝 . We describe the workflow below. For simplicity we let 𝜔 = [0, 𝐿] 3 . Let 𝑁 be a fixed positive integer. We partition 𝜔 into 𝑁 3 hexahedral cells ω𝑛 , (𝑛 = 1, 2, … , 𝑁 3 ), numbered with 𝑛 rather than 𝑖𝑗𝑘 for brevity, and with uniform side length ℎ = 𝐿∕𝑁. (𝑓 ) and 𝜔 (𝑟) We now describe which cells should be considered part of 𝜔 (𝑓 ,′) or 𝜔 (𝑟,′) . For the rules of this assignment, there are several factors to consider. First, the discrete porosity |𝜔 (𝑓 ,′)
|∕|𝜔| should be close to the true porosity |𝜔 (𝑓 ) |∕|𝜔|. Second, the connectivity of the discrete pore space should be similar to the connectivity of the true pore space. This second point becomes critical for assemblies in 𝑑 = 3, as the complexity of the pore-space is substantial. Finally, the assignment rule should be consistent as we refine the grid, i.e., make ℎ smaller. (𝑓 ) and 𝜔 (𝑟) described in Section 3.1. Gray cells represent 𝜔 (𝑟,′) and white cells represent 𝜔 (𝑓 ,′) . Contours of three particles 𝐵 𝑝 are shown in black. To compare, it is useful to study in each scenario the assignment of the cell ω which is third from left in the top row.
We now describe the assignment rules, i.e., how 𝜔 (𝑟,′) and its complement 𝜔 (𝑓 ,′) are defined. We illustrate several scenarios of assignment (I, II, . . . , VI) in Fig. 1 for a coarse grid, and in Fig. 2 for a fine grid. Formally, we let
This rule corresponds to rule I in both figures, with the complementary definition
After testing the strategies illustrated in Fig. 2, we found the strategy I to be the most robust and versatile. In particular, strategies II-III overpredicted while IV underpredicted the porosity, and V-VI lacked symmetry.
However, after this first step of the assignment of cells, there might be so-called ''dead'' pores that are not connected to other pores, and do not play a role in the global flow solution. This is handled next. (𝑟,′′) and 𝜔 (𝑓 ,′′) by trimming and padding
The construction of 𝜔 (𝑓 ,′) described above creates generous void space around the space occupied by the particles, since we start with a bounding box designed to contain all the particles, and they are not always arranged close enough to the sides. Solving for the flow in a domain with large void space around the particles would direct the flow around this domain and would give artificially large conductivity after upscaling. To prevent this artifact, we trim a small percentage from each face of 𝜔. Let 𝑝 1 ∈ [0, 1) such that 𝑝 1 𝑁 is an integer (our choice was 𝑝 1 = 0.02) and set
We then set
In flow solutions, it is common to choose a parabolic inflow boundary condition. Following our prior experience in [22] we apply the following strategy to make this assignment easy as well as to make the flow solutions comparable between different synthetic geometries. We pad the trimmed 𝜔 𝑡𝑟𝑖𝑚 by void space to use the portion of the new void space as the inlet boundary. Let 𝑝 2 ≥ 0, such that 𝑝 2 𝑁 is an integer number of cells we wish to pad with (our choice was 𝑝 2 = 0.05) and set
Z. Hilliard et al.
Follow-up on illustration in Fig. 1: working with a fine grid to assign the cells to 𝜔 (𝑓 ) and 𝜔 (𝑟) . Gray cells represent 𝜔 (𝑟,′) and white cells represent 𝜔 (𝑓 ,′) .
Then 𝜔 𝑝𝑎𝑑 can be discretized into 𝑀 = 𝑝 2 (1 -2𝑝 1 ) 2 𝑁 3 cells ω𝑛 (𝑛 = 𝑁 3 + 1, … , 𝑁 3 + 𝑀) with sidelength ℎ. We set
It is possible that 𝜔 (𝑓 ,′′) found in Section 3.2 is not connected. This would make solving the Stokes problem in 𝜔 (𝑓 ,′′) ill-posed. To circumvent this problem, we take another step in approximation of geometry in which we remove the cells not connected to the inflow and outflow boundaries of 𝜔 (𝑓 ,′′) .
For example, define the inflow and outflow boundaries
We will now explain how to modify 𝜔 (𝑓 ,′′) by introducing the notion of orthogonally connected cells. We will say that a cell ω𝑖 ⊂ 𝜔 (𝑓 ,′′) is orthogonally connected to ω𝑗 ⊂ 𝜔 (𝑓 ,′′) if there is a sequence of cells ω𝑖 = ω𝑛 1 , ω𝑛 2 , … , ω𝑛 𝑚 = ω𝑗 contained in 𝜔 (𝑓 ,′′) such that each pair ω𝑛 𝑙 and ω𝑛 𝑙+1 shares a common face for 𝑙 = 1, 2, … , 𝑚 -1. We say that a cell ω𝑖 is orthogonally connected to 𝛤 𝑖𝑛 if there is at least one cell ω𝑗 ⊂ 𝜔 (𝑓 ,′′) such that one face of ω𝑗 is in 𝛤 𝑖𝑛 and ω𝑖 and ω𝑗 are orthogonally connected. We now define ′′) and ω𝑛 is orthogonally connected to
An illustration of the concerns on connectivity is given in Fig. 3. The hatched cells are those in 𝜔 (𝑓 ,′′) but not in 𝜔 (𝑓 ) as they are not orthogonally connected to 𝛤 𝑖𝑛 .
Now we illustrate the first two steps in (5), first in a cartoon in Fig. 3. Next, Fig. 4 shows the set 𝜔 (𝑟) = ⋃ 𝑝 𝐵 𝑝 , as well as its discretization into hexahedra, followed by the result of trimming and padding the domain, and finally by removal of any dead cells to obtain 𝜔 (𝑟,′′) .
In this Section we describe how to conduct the flow simulations in step (5c) following (5a)-(5b). (𝑓 ) in white, rock cells in 𝜔 (𝑟) in gray; the hatched cells are the fluid cells that must be eliminated from 𝜔 (𝑓 ) since they are not connected to 𝛤 𝑖𝑛 . Bottom: illustration of the process of obtaining 𝜔 (𝑓 ) in 2D. From left: 𝜔 (𝑓 ,′) found using strategy I from Fig. 1, 𝜔 (𝑓 ,′′) obtained by trimming and padding, and the final domain 𝜔 (𝑓 ) obtained by elimination of dead-end cells. Fig. 4. Illustration of first two steps in (5). From left to right: synthetic geometry as described in Section 2, discretization into hexahedra (𝜔 (𝑟,′) ) as described in Section 3.1, and 𝜔 (𝑓 ,′′) found after trimming and padding the discretization (𝜔 (𝑟,′′) ) as described in Section 3.2.
To model the fluid flow through the pore space 𝜔 (𝑓 ) we use Stokes flow model given as
where 𝑈 is the pore-scale velocity, 𝜇 𝑓 is viscosity, 𝑃 is the pore-scale pressure, 𝜈 is the outward normal unit vector, 𝜔 (𝑓 ) is the fluid domain discretized into hexahedral cells.
The flow model (16a)-(16b) is supplemented with boundary conditions (16c)-(16e). It remains to identify the inflow, outflow, and no-flow boundaries 𝛤 𝑖𝑛 , 𝛤 𝑜𝑢𝑡 , 𝛤 0 , respectively.
In this paper we conduct simulations of ( 16) for the needs of upscaling and thus require flow solutions corresponding to at least three distinct flow solutions associated with each 𝜔 (𝑓 ) . We choose these flow directions to approximately align with 𝑥 1 or 𝑥 2 , or 𝑥 3 directions. For each of these we set up the boundary conditions and simulate the flow, and thus obtain the solutions 𝑈 (𝑖) , 𝑃 (𝑖) for 𝑖 = 1, 2, 3.
To be specific let 𝑖 = 1. We set 𝛤 𝑖𝑛 as the inflow boundary that corresponds to the face of 𝜔 (𝑓 ) with the least 𝑥 1 value, 𝑈 𝐷 = (𝑢 𝐷 , 0, 0) is a specified parabolic profile in the 𝑥 1 -direction, 𝛤 𝑜𝑢𝑡 is the outflow boundary that corresponds to the face of 𝜔 (𝑓 ) with greatest 𝑥 1 value, and 𝛤 0 = 𝜕𝜔 (𝑓 ) ⧵ 𝛤 𝑖𝑛 ⧵ 𝛤 𝑜𝑢𝑡 . An illustration in 𝑑 = 2 of these sets is given in Fig. 3.
We follow an analogous process in the directions 𝑖 = 2 and 𝑖 = 3.
The hexahedral discretization of 𝜔 (𝑓 ) lends itself to the well known marker and cell (MAC) scheme which delivers approximations to 𝑈 and 𝑃 . A detailed description of the MAC scheme is given in [51], with a recent Ref. [52] applying MAC scheme to a more complex problem than Stokes flow. It is also worth noting that the MAC scheme produces conservative approximations to 𝑈 which are in the space 𝑅𝑇 [0] popular in mixed finite element setting on hexahedral grids [53], with the approximations to 𝑃 in the space of piecewise constants over the grid ℎ (𝜔 (𝑓 ) ).
The MAC scheme introduces the degrees of freedom (DOFs) on the faces and centroids of each ω𝑛 ⊂ 𝜔 (𝑓 ) in the following manner. The DOFs for 𝑃 are located at the centroid of each ω𝑛 and the DOFs for 𝑈 𝑖 (the 𝑖th component of 𝑈 ) are located on the faces of each ω𝑛 that are orthogonal to the 𝑖th coordinate direction. After the finite difference approximation is written out for each component of 𝑈 on a staggered grid involving grid edges, we obtain a linear system which is solved by an iterative solver with state-of-the-art preconditioners [23]. The numerical solution (𝑈 ℎ ) is then extended to 𝜔 (𝑟,′′) and the rest of 𝜔 (𝑓 ,′′) by 0 on any cell not contained in 𝜔 (𝑓 ) .
In this paper we use the implementation of the MAC scheme adapted to the pore-scale geometry described in [23] and implemented in HybGe-Flow3d [54].
We illustrate the solution to (16) in Fig. 5. This solution is next processed for upscaling as discussed in Section 5.
In this Section we describe how to conduct the upscaling step (5d) following (5a)-(5c). We assume (5c) is completed, and the numerical approximations to (16) are available. We post-process these to obtain the conductivity 𝐾 = 𝜇 -1 𝑓 𝜅, and give results. Assume we have the approximations (𝑈 ℎ,(𝑖) , 𝑃 ℎ,(𝑖) ), 𝑖 = 1, 2, 3 to the Stokes problem (16) as described in Section 4, with 𝑖 = 1, 2, 3 corresponding to each of the coordinate directions. We use their averages over the REV to calculate 𝜅. Below we describe the pre-processing of (𝑈 ℎ,(𝑖) , 𝑃 ℎ,(𝑖) ), 𝑖 = 1, 2, 3, next the actual method of averaging and fitting the results to the conductivity, followed by the notes on scaling of the resulting conductivities. In Section 5.4 we present the results. 5.1. Restricting (𝑈 ℎ,(𝑖) , 𝑃 ℎ,(𝑖) ), 𝑖 = 1, 2, 3 to avoid numerical artifacts Here we aim to remove various numerical artifacts that might be present in 𝑈 ℎ, (𝑖) . In particular, the numerical solutions near the boundary of 𝜔 𝑡𝑟𝑖𝑚 may exhibit certain boundary effects such as circulating flow near 𝛤 𝑜𝑢𝑡 or the impact of setting 𝑈 = 0 on the lateral boundaries of 𝜔 (𝑓 ) . In order to eliminate as many of these effects as possible, we restrict the solutions (𝑈 ℎ,(𝑖) , 𝑃 ℎ,(𝑖) ) defined on 𝜔 (𝑓 ) to
Here 𝑝 3 ∈ [0, 1) is the fraction to further inset the domain; in practice, we set 𝑝 3 = 0.09. This strategy was first applied in [22] and refined in [24] where we found that it helps to avoid numerical artifacts. Regarding the modeling accuracy of this strategy, in [22,24] we found that the conductivities 𝜅 depend mildly on the size of REV, i.e., on the value of 𝑝 3 . The choice 𝑝 3 = 0.09 is a compromise between the need to avoid numerical artifacts and the aim to keep the REV as large as possible.
First, recall Darcy's law (1c) with ∇𝐷 = 0,
where 𝐾 is conductivity and 𝑣 and 𝑝 are the Darcy scale velocity and pressure respectively. This law can be derived via homogenization of the Stokes equation [18][19][20][21]. The process of homogenization involves periodic boundary conditions and an impulse forcing term 𝑒 𝑖 (the 𝑖th standard basis vector in R 𝑑 ) in (16a) as well as averaging of ∇𝑈 . However, quality results can be obtained with a simpler process for the flow driven by inflow boundary conditions (i.e., (16) as written) [22,55]. We follow the upscaling method from [55] which we now summarize. Given the solutions (𝑈 ℎ,(𝑖) , 𝑃 ℎ,(𝑖) ) to ( 16) over 𝜔 𝑖𝑛𝑠𝑒𝑡 , we identify their averages with 𝑣 and ∇𝑝 as follows. Fix the coordinate direction 𝑖, and bisect 𝜔 𝑖𝑛𝑠𝑒𝑡 perpendicularly to 𝑥 𝑖 into sets 𝜔 𝐿 𝑖 and 𝜔 𝑅 𝑖 such that:
(1) Each 𝜔 𝐿 𝑖 and 𝜔 𝑅 𝑖 is a rectangular prism.
(2) For each 𝑖, 𝜔 𝐿 𝑖 and 𝜔 𝑅 𝑖 are disjoint with the exception of sharing a common face.
(3) For each 𝑖, 𝜔 𝐿 𝑖 ∪ 𝜔 𝑅 𝑖 = 𝜔 𝑖𝑛𝑠𝑒𝑡 (see Fig. 6).
We now define 𝛿𝑋 𝑖 as the distance in the 𝑖th coordinate direction between the centroids of 𝜔 𝐿 𝑖 and 𝜔 𝑅 𝑖 and define the averages
In general each of the terms in ( 19) is a vector of 𝑑 components. Now we use ( 18) to obtain the system of 𝑑 × 𝑑 linear equations
which can then be solved for 𝐾 𝑗𝑚 . Here 𝑣 (𝑖) 𝑗 is the 𝑗th component of 𝑣 (𝑖) and 𝐾 = (𝐾 𝑗𝑚 ). For the experiments we described, typically, 𝑣 (𝑖) 𝑗 as well as ∇ ℎ 𝑝 (𝑖) 𝑗 are close to zero when 𝑖 ≠ 𝑗. Thus, while the full tensor 𝐾 could be computed, in this paper we only calculate and report only the main diagonal values 𝐾 11 , 𝐾 22 , 𝐾 33 of 𝐾 found from (20) for 𝑖 = 𝑗 = 1, 2, 3. This approach is sufficient for geometries that are either nearly isotropic or for which the anisotropy is aligned with the coordinate axes.
We make several observations which help to set up simulations of a physical system in non-dimensional units and translate these later to a physical set-up with a fixed grain size so that we can obtain a corresponding permeability value given in the usual units [m 2 ]. In our prior experience from [22,24] the process delivers the values 𝜅 agreeing with typical permeability values found experimentally for porous media with the given grain size. We note that the methods outlined in Sections 4 and 5.2 are most useful to find the relative conductivity in response to deformation; these relative values are independent of the particular scaling process.
Let (𝑈 , 𝑃 ) be the solution to ( 16) and ( Ũ , P ) the solution when the boundary condition (16c) is replaced with Ũ𝐷 = 𝑐𝑈 𝐷 for some 𝑐 ≠ 0. One can verify that Ũ = 𝑐𝑈 and P = 𝑐𝑃 . If the solution ( Ũ , P ) is upscaled by the method in Section 5, then one arrives at the system of Eqs. (20) with both sides multiplied by 𝑐 so that the same upscaled conductivity K = 𝐾 is found.
Let (𝑈 , 𝑃 ) be the solution to ( 16) and ( Ũ , P ) the solution to the flow model -μ 𝑓 𝛥 Ũ -∇ P = 0 replacing (16a), with the boundary condition Ũ𝐷 = (𝜇 𝑓 ∕μ 𝑓 )𝑈 𝐷 in (16c). One can verify that Ũ = (𝜇 𝑓 ∕μ 𝑓 )𝑈 and P = 𝑃 . Upscaling the solution ( Ũ , P ) results in the conductivity K = (𝜇 𝑓 ∕μ 𝑓 )𝐾.
Let (𝑈 , 𝑃 ) be the solution to ( 16) and ( Ũ , P ) the solution on ω(𝑓) which is congruent to 𝜔 (𝑓 ) but scaled by a factor 𝑐 = L∕𝐿 for some L > 0 (i.e., Ũ and P are defined on x ∈ ω(𝑓) = 𝑐𝜔 (𝑓 ) and we use the boundary condition Ũ𝐷 (x) = 𝑈 𝐷 (x∕𝑐), which corresponds to a change of variables x = 𝑐𝑥). Then Ũ (x) = 𝑈 (x∕𝑐) and P (x) = 𝑃 (x∕𝑐)∕𝑐. Upscaling the solution ( Ũ , P ) produces the conductivity K = 𝑐 2 𝐾 after noting that 𝛿 X𝑚 = 𝑐𝛿𝑋 𝑚 so that ∇ ℎ p(𝑖) 𝑚 = ∇ ℎ 𝑝 (𝑖) ∕𝑐 2 .
We apply the process described above in Sections 5.1-5.3 to the assemblies found in Section 3 for each of the scenarios in Section 2.2. We obtain conductivities for each of these scenarios, with focus on the relative changes in the conductivities associated with the scenarios of the geometries undergoing a loading and unloading cycle.
As a by-product of pre-processing step which delivers 𝜔 (𝑓 ,′′) described in Section 3 we obtain the porosities 𝜂. The changes in 𝜂 for the 𝑘 = 1, … , 𝐾 scenarios can be correlated to the relative conductivities.
In upscaling, we set the length scale to 𝐿 = 1 [cm so that the particles are on the order of 10 -4 [m] in diameter and the computed upscaled conductivities are in the units of [Pa -1 s -1 cm 2 ].
Below we present an overview followed by a detailed discussion; the raw data tables are given in the Appendix in Appendix A.1. In Section 5.4.3 we conclude with a discussion how the upscaled data can be used in a Darcy scale model.
An overview of results is given in Fig. 7, with detailed results given in Section 5.4.2. First, we see that the process we described in (5) with details in Sections 2-5.3 is robust. We also see, as expected, that the conductivities have order of magnitude reasonably matching the conductivities measured experimentally [17]. Furthermore, generally, conductivities and porosities decrease with load by about maximum 20% within the range considered. Additionally, the conductivity tensors corresponding to these geometries are essentially isotropic (insofar as it is difficult to definitively say that they are anisotropic). These observations are consistent across all the scenarios in Section 2.2.
The values across the different scenarios generally match, with the exception of scenario (C) which features conductivities smaller by a factor of about 2.2. Recalling that this scenario involves more particles, the smaller conductivity can be explained by the fact that we are packing more particles to the roughly similar box size. This is exacerbated by the use of the same fix grid resolution 𝑁; on the other hand, scenario (C) generally enables tighter packing, somewhat smaller porosity, and more tortuous flow paths which result in lower conductivity.
Lastly, we discuss a certain common trend present in all results. For example, consider results plotted in Fig. 13 for scenario (D). Here we can see hysteresis in the conductivity, i.e., values following a different path at loading than at unloading. It appears however that this hysteresis is arises in the load-porosity relationship rather than appears as standalone hysteresis in the load-conductivity relationship. The phase plane of porosity vs conductivity shows much less sensitivity to the direction (loading or unloading). 6. Markers in white denote unloading.
Now we provide detailed discussion of the upscaling results.
Results from the experiment (A) described in Section 2.2.1 are plotted in Fig. 14 and given in Table 5. Overall, as expected, we see that the conductivities and the porosity decrease with load. In addition, we see hysteresis in the values when we switch from loading to unloading. This observation of hysteresis is consistent across all experiments.
Z. Hilliard et al. 9. Markers in white denote unloading.
In Fig. 8 we give plots for (B1) described in Section 2.2.2, with data in Table 6. We see the same hysteretic trend as (A), but it is more exaggerated. Additionally, there is an outlying data point for 𝐾 22 at 500 [kPa] on the unloading cycle. We hypothesize that there is some local irregularity in the geometry that caused difficulty for the Stokes solver for flow in the 𝑥 2 -coordinate direction.
The results for (B2) described in Section 2.2.3 are plotted in Fig. 9, with data in Table 7. The conductivity values in this experiment are more tightly grouped than in (B1). However, the Stokes solver failed to converge for the geometry at 800 [kPa] on the loading cycle for all flow orientations. 10. Data from the geometries in experiment (D) that were used to generate the 𝛱 (𝑖) geometries is included for comparison. In Table 8 we see the results for (C) described in Section 2.2.4, in which we increased 𝑁 𝑝 to 24,000. Increasing the number of spheres decreased conductivities and porosities as compared to other experiments, but we see the same qualitative trends. On the other hand, the iterative solver for the Stokes flow model failed to converge for the geometry at 200 [kPa] on the unloading cycle for all flow orientations. Results from this experiment are plotted in Fig. 10. Lastly, we discuss the experiment (D) from Section 2.2.5. The general trends that are seen in the other experiments are also seen here. In particular, there is an interesting ''bump'' in the conductivities in the 700-800 [kPa] range for both the loading and unloading cycles. Results from this experiment are plotted in Fig. 11.
We continue next with data from (D), but we remove some of the particles; the case is recorded as (D.r) described in Section 2.2.6. Removing spheres results in an increase of the porosity as well as of the conductivity as one should expect; see the results plotted in Fig. 12. 5. Markers in white denote unloading.
Fit of conductivity to 𝐾 = 𝐶 1 𝜂 𝛼 1 (non-normalized) 𝐾 = 𝐶 2 (𝜂∕𝜂 0 ) 𝛼 2 (normalized). Normalization was done to data at load 100 [kPa] on the loading (L) cycle. The data for the (D.r) fit was normalized by the same data as the (D) case. The cases are described in Section 2.2 and fitted to the data discussed in Section 5.4.2.
𝐶
In the end, we need the upscaled conductivities to be useful in a Darcy model for flow and deformation. In particular, it is not practical to imagine we could be solving a new (5a)-(5d) experiment each time a load changes at the Darcy scale so as to determine the new conductivities.
For this reason, we collect the results in a look-up table such as the Tables in Appendix A.1. Furthermore, we attempt a fit of the data in these look-up tables.
We focus here on the porosity-conductivity relationship which follows a power law (i.e., Carman-Kozeny) for each geometry, as can be seen in Fig. 15 and Table 2. It is interesting to note that the exponent in Carman-Kozeny relationship (4) does not change much when the spheres are removed as described in Section 2.2.6 and shown in Fig. 12 and Fig. 15.
Thus the load-conductivity data can be used as follows: we determine porosity for a new load, then find the conductivity for that porosity from Carman-Kozeny law or from a look-up table.
In this section we illustrate how to use the data found in the steps (5a)-(5d) of ( 5) in a Darcy scale model for flow and deformation.
To this aim we consider the linear Biot model of poro-elasticity which we modify. In particular, we introduce the dependence of the conductivity 𝐾 = 𝜇 -1 𝑓 𝜅 on the deformation, and in fact we account for this through the porosity 𝜂 as an intermediate step. The quasi-static Biot model is well known. See [56](Pg. 114) for the formulation and [57] for the well-posedness. For notation, we follow closely [58]. We recall now the model (1) for ease of exposition neglecting the gravity terms.
We let 𝛺 = (0, 1) ⊂ R represent a porous medium that is fully saturated with water. We recall the model (1) for convenience, and we set ∇𝐷 = 0, 𝑓 = 0 in (1) focusing on the influence of 𝑓 𝑓 ≠ 0. The elastic deformation and flow equations read, respectively,
where 𝑢 is the displacement
where 𝜅 is the permeability [m 2 ] of the media and 𝜇 𝑓 is the viscosity [Pa s] of liquid water. The physical parameters used in ( 21) are given in Table 1. We also have the Biot constant 𝛼 which as usual we set 𝛼 = 1 [-]. We shall also impose the following boundary conditions
Finally, the model ( 21) requires also initial conditions on the fluid content 𝑐 0 𝑝+𝛼∇⋅𝑢 which we simplify assuming initial equilibrium ∇ ⋅ 𝑢| 𝑡=0 = 0 and require
Under the assumptions that 𝜆, 𝜇, 𝜇 𝑓 , 𝑐 0 are positive, and that 𝜅 is a symmetric uniformly positive definite tensor, the problem ( 21) is well posed [57].
We focus now on the dependence of 𝜅 on the mechanical deformation. First, we replace the direct dependence of 𝜅 on the deformation, e.g., on strain 𝜅 = 𝜅(∇ ⋅ 𝑢) by
Table 3 Results of Example 6.1 found with the numerical model (25a) and 𝑀 = 50, 𝜏 = 0.01. We present the magnitude of displacements and pressures, impact on porosity at 𝑡 = 3 when the injection stops. Also, we present the time 𝑇 𝑒𝑞 to equilibrium defined as the earliest time 𝑡 𝑛 when max 𝑥 𝑝(𝑥, 𝑡 𝑛 ) < 1. Generally the more dramatic modification corresponds to the more influence on 𝑝. The impact on the displacements 𝑢 and the porosities 𝜂 is milder but also visible. For (22) to be useful, we need 𝜂(𝑥, 𝑡) as well as the actual relationship 𝜅(𝜂) itself. The latter is identified as part of results in Section 5.4.3.
The former, the knowledge of the current porosity 𝜂(𝑥, 𝑡) follows from the porosity equation given in [59](Eq. ( 18)), [7](eq.( 3b)).
Remark 6.1. We note that the porosity value given in ( 23) is the macroscopic variable relevant at the Darcy scale; this quantity changes locally with the deformation. The porosity variable maintains the interpretation as the local value |𝜔 (𝑓 )
|𝜔| but this connection to any specific 𝜔 (𝑓 ) (𝑥) or to the local deformation of particles found by DEM is not maintained.
Numerical schemes and simulations for the Biot model have been the subject of intense research for a variety of applications. See e.g. the results in [2,4,5,12] for some of the numerical schemes and simulation results which have directly informed our work in this paper.
To define the discrete Biot model, we use the well known three-field formulation which approximates 𝑝, 𝑞 𝑓 using mixed finite element method of lowest order 𝑅𝑇 [0] implemented as cell centered finite differences (CCFD). For the displacements we use piecewise linears. These choices are the same as in [2,4]. For implementation we adapt the code in our Biot capsule [60].
Let 𝑈 𝑛 , 𝑃 𝑛 , 𝑄 𝑛 𝑓 the degrees of freedom of the approximations to 𝑢, 𝑝, 𝑞 𝑓 at time 𝑡 𝑛 to the solutions to (21). We seek these as the solutions to
where the matrices {𝐴 𝑢𝑢 , 𝐴 𝑝𝑢 , 𝐴 𝑞 𝑓 𝑝 , 𝑀 𝑝𝑝 , 𝑀 𝑞 𝑓 𝑞 𝑓 } and vectors { 𝑛 , 𝑛 , 𝑛 } are found in the usual way; see, e.g., [58].
In particular, the matrix 𝑀 -1
is made of the cell-wise transmissibilities directly computed from the conductivities 𝜅 𝜇 𝑓 arranged on a diagonal. Thus the entries of matrix 𝑀 -1
𝑞 𝑓 𝑞 𝑓 are proportional to 𝜅, i.e., in turn, depends on 𝑈 . We denote this by writing 𝑀 𝑞 𝑓 𝑞 𝑓 (𝑈 𝑛 ) in (24). In the end the problem ( 24) is nonlinear.
To close, we now require the particular relationship (3) as well as a solver capable of solving the resulting system (24). As we said we will follow a sequential approach as well as (22).
In this paper we take the following approach. We resolve the nonlinearity in (24) by time-lagging in a sequential approach, with the porosity treated as an auxiliary variable. We linearize (24) by time-lagging and replacing (𝑈 𝑛 ) by (𝜂 𝑛-1 ).
Concretely, we proceed as follows. When entering time step 𝑛, we assume 𝜂| 𝑡 𝑛-1 is known. (When 𝑛 = 0 we use 𝜂 0 ). With this we evaluate 𝜅 𝑛-1 = 𝜅(𝜂 𝑛-1 ) cell-wise. From this step we obtain the transmissibilities and recalculate the 𝑀 𝑞 𝑓 𝑞 𝑓 (𝜅 𝑛-1 ) component of (𝜂 𝑛-1 ). Next we solve the linearized version of ( 24) Once the new values (𝑃 𝑛 , 𝑈 𝑛 , 𝑄 𝑛 𝑔 ) are known, we update the value of porosity from the time-discrete explicit approximation to ( 23)
which is straightforward, since the displacements in this 𝑑 = 1 model have degrees of freedom at the grid's vertices, and porosities are associated with the cell centers.
With the new cell-wise values of 𝜂 𝑛 known, we evaluate
to be used in the next time step. Here we use a look up table or a Carman-Kozeny relationship that we obtain in Section 5.4.3. 4) and the exponent 𝛼 𝐶𝐾 = 0, 5, 7 (from top to bottom). We see a significant impact in the conductivity 𝜅 and some impact on the pressure 𝑝 and displacement 𝑢.
In principle, one could iterate on 𝜂 in (25a) to obtain a consistent solution to (24). Alternatively, one could also implement a faster solver such as Newton's method. This approach is however very complicated, and not warranted by the typical characteristic time of flow (fast) versus the deformation (slow).
We illustrate the strategy given above by a suite of examples. We do not aim to use any particular DEM results, but rather show how the simulations are affected by the choice of data, and in particular, by the parameters 𝛼 𝐶𝐾 in (4). ). We impose 𝑢| 𝜕𝛺 = 0, and 𝑝(0, 𝑡) = 0, 𝜅∇𝑝 ⋅ 𝜈| 𝑥=1 = 0, and start from initial equilibrium with 𝑝 𝑖𝑛𝑖𝑡 (𝑥) = 0, 𝑢 𝑖𝑛𝑖𝑡 = 0. The dynamics is driven by fluid injection from 𝑓 (𝑥, 𝑡) = 10 for 𝑡 ∈ (0, 3]. For 𝑡 > 3 there are no sources, but the fluid can drain freely through the boundary at 𝑥 = 0.
For numerical discretization we employ a uniform grid with ℎ = 1 𝑀 , and uniform time-stepping with 𝜏 = 0.1 or less. We first run the simulation to understand the dynamics, and to test the grid convergence.
As expected, we find that during the initial transients 𝑡 < 0.3 we see an increase of the pressure 𝑝(𝑥, 𝑡) close to the injection region, and the displacements and porosities react appropriately. This means that the injected fluid will cause pressure to increase and displacements to vary until equilibrium is attained around 𝑡 = 3. Then the transients stabilize, and about 𝑡 = 3 the solutions are close to some steady state: the injected fluid is almost completely drained at any time through the boundary 𝑥 = 0. We turn off the fluid source at 𝑡 = 3, and the fluid now drains until the pressures and displacements stabilize towards the new steady state. At 𝑡 > 9 the pressures are very close to 𝑝 ≈ 0. In these simulations we do not allow 𝜅 to change.
We also study sensitivity of the model to the choice of spatial and time grid. We choose 𝑀 = 50,100 and various time steps 𝜏. In Fig. 16 we present the results of 𝑀 = 100, 𝜏 = 0.001 at three separate time steps 𝑡 = 0.2, 𝑡 = 3, 𝑡 = 5 corresponding to initial transients, quasi-equilibrium, and long-term equilibrium. The results show that for a grid finer than 𝑀 = 100, 𝑡 = 0.001 the numerical solution is almost the same as that for 𝑀 = 100, 𝑡 = 0.001. In other studies we choose therefore a slightly coarser grid 𝑀 = 50, 𝑡 = 0.01. 4) with exponents 𝛼 𝐶𝐾 = 5 and 𝛼 𝐶𝐾 = 7 and we compare these to that when 𝜅 is fixed. Results are given in Fig. 17 and in Table 3.
Based on this example we see some dependence of the solutions {𝑢, 𝑝, 𝜅} on the choice of the model for the altered permeability 𝜅 = 𝜅(𝜂), with the largest impact on the pressures. We see that the pressure values are sensitive to the specific modification of conductivity 𝜅 with difference of up to 5%, while the displacements and porosities differ by about 1% to 0.1%, respectively. We also see some impact on the time to equilibrium included in Table 3. Example 6.3. Finally we address the issue of stiffness, i.e., the sensitivity of the model to the assumed Lame parameters 𝜆, 𝜇. We alter these in the simulation setup from Example 6.1; recall these are 𝜆 0 = 2.8846 × 10 7 [Pa] 𝜇 0 = 1.9231 × 10 7 [Pa]. Then we set 𝜆 1 = 1.5 × 10 7 , 𝜇 1 = 10 7 [Pa] corresponding to about twice smaller 𝐸, solve the problem for the corresponding, 𝑢 (1) , 𝑝 (1) as well as porosity and permeability changes. Next we consider a ''stiffer'' case where 𝜆 2 = 610 7 , 𝜇 2 = 410 7 [Pa] are about twice stiffer than 𝜆 0 , 𝜇 0 , respectively, and solve the problem for 𝑢 (2) , 𝑝 (2) . The results are given in Table 4.
The results of Example 6.3 are qualitatively similar to those in Example 6.1, thus we do not plot these. However, Table 4 shows that there is a noticeable quantitative difference between the displacements, pressures, and porosity change. In particular, as expected, the change of stiffness by a factor of 2 results in an inversely proportional dependence of the displacement, with an interesting significant difference on the porosity changes not correlated to the factor 2 or 0.5. In summary, we see that the Biot model responds to the altered permeability as well as to the changes in mechanical parameters. In general, the pressures respond more strongly to the permeability changes, and the displacements respond more strongly to the stiffness changes. More studies are needed in the future.
In this paper we presented our efforts on connecting the pore-scale models of flow and deformation to the Darcy scale. Specifically, we considered generation of physically realistic pore-scale geometries in response to some mechanical load using DEM. Next we applied a pore-scale flow model and obtained upscaled values of conductivity 𝜅 corresponding to these loads. In the end we used this data in a Biot model which has been altered to allow the incorporation of the evolving 𝜅.
We have various conclusions. First, we believe that our process and first attempt to connect DEM, upscaling and Biot model, was successful. However, it also came with various challenges. We simulated physical realistic assembles of particles with DEM, but, while we wanted to include as many particles as possible, we faced restrictions from subsequent steps. Now, as in any upscaling, we aim for the size of REV 𝜔 to be adequate and large. However, we are limited by the resolution of flow simulations. In particular, to resolve the geometry around the particles appropriately, we required a grid of 200 3 over 𝜔; such resolution creates difficulties for the flow solver, requiring a robust iterative large scale linear solver able to handle about #DOF ≈ 𝑂(𝜂 × 10 7 ) unknowns. Since we are not able at present to simulate Stokes flow on significantly larger grids than with this #DOF, these constraints restrict the number 𝑁 𝑝 of particles in the assemblies to 𝑁 𝑝 = 𝑂 (10 4 ). On the other hand, this 𝑁 𝑝 seemed sufficient from the mechanical standpoint and produced reasonably dependable predictions.
Second, the originally linear Biot model becomes nonlinear when we use the relationship (3) between conductivity and deformation. At the same time the DEM simulations describe the evolution of the pore-scale due to the deformation; the DEM results we see (and use) do not quite fit within an elastic response unlike that prescribed by the Biot model. In particular, we see hysteretic behavior apparent in the porosity and conductivity changes due to the load which cannot be accounted for by our current version of the Biot model. We also note an apparent sensitivity of the solutions to the mechanical properties assumed in the simulations. These features suggest the need to define a process to obtain more constitutive data from DEM or other methods consistent with a particular DEM generated assembly of particles and the considered scenarios.
We plan to address the challenges and questions raised above in future work.