FLAC3D Theory and Background • Fluid-Mechanical Interaction

One-Dimensional Filling of a Porous Region

Note

To view this project in FLAC3D, use the menu command Help ? Examples…. Choose “Fluid/1DFillingPorousRegion” and select “1DFillingPorousRegion.f3dprj” to load. The project’s main data files are shown at the end of this example.

In this problem, flow is driven through an initially dry porous layer of large lateral extent under a constant pressure, \(p_{_{o}}\), applied at the base. The transient location of the filling front is compared to an exact sharp-front solution for the cases with and without gravity.

Voller et al. (1996) give an analytic solution for this problem under the assumptions of a sharp-front, rigid-porous matrix and incompressible Newtonian fluid. In their solution, the flow is governed by Darcy’s law, and there is a constant atmospheric pressure in the air ahead of the free surface.

Let the \(x\)-axis of reference be oriented in the direction of flow, with the origin at the base of the layer. The solution for the front location, \(x_{f}\), may be expressed in terms of two dimensionless variables: \(\widehat{t}=\) \(t/T^{\ast }\)and \(\widehat{x}=x_{f}/L^{\ast }\); and a dimensionless parameter, \(\gamma\), defined using the expressions \(T^{\ast}=n\mu /p_{o}\), \(\ L^{\ast }=\sqrt{\kappa }\) and \(\gamma =\sqrt{\kappa }\) \(\rho g/p_{o}\). In these equations, \(n\) is porosity, \(\kappa\) is intrinsic permeability (product of mobility coefficient, \(k\), and dynamic viscosity, \(\mu\)), \(\rho\) is fluid density, and \(g\) is gravity.

When gravity is ignored, the solution has the form

(1)\[\widehat{t}={{1}\over{2}} \left[ \widehat{x} \right] ^{2}\]

For filling under gravity, the front location is given by

(2)\[\widehat{t}=-{\widehat{x}\over\gamma} -{\ln \left( 1-\gamma \widehat{x}\right) \over\gamma ^{2}}\]

Equation (2) may be shown to converge to the no-gravity solution when \(\gamma\) goes to zero.

The numerical solution to the filling problem is presented in dimensionless form. To derive this solution, scaled properties are used in the simulations: \(p_{o}\) = 1, \(k\) = 0.25, \(n\) = 0.5; and for gravity flow, \(\rho _{w}\) = 1 and \(g\) = 1. Further, using \(\mu\) = 4 in the preceding definitions, the characteristic parameters for the simulation are \(T^{\ast }\) = 2, \(L^{\ast }\) = 1, \(\gamma\) = 1.

To simulate incompressible flow, the bulk modulus of the fluid, \(K_{f}\), is assigned a value that is large compared to the pore-pressure variations in the simulation (\(K_{f}\) = 100). For both cases, the grid corresponds to a column of 25 zones, 0.625 units high and 0.025 units wide. The initial value of pore pressure and saturation is zero. The pore pressure is fixed at \(p_{o}\), and the saturation is given a value of 1 at the base of the model. The simulation is conducted for a total of 0.25 units of time (\(\widehat{t}\) = 0.125). A FISH function, flacfront, captures the times at which nodes reach a saturation of 1% and 99%. The analytic sharp-front solution is evaluated by another FISH function, solution. The data file for the case without gravity is listed in “1DFillingPorousRegion-NoGravity.f3dat”; the data file for the case with gravity is listed in “1DFillingPorousRegion-Gravity.f3dat”.

The results with and without gravity are presented in Figure 1 and Figure 2, respectively. As seen in these figures, the sharp-front solution is bounded above and below by the 99% and 1% saturation fronts. In fact, the vertical distance between these fronts corresponds directly to the grid size in the direction of propagation of the filling front. (The saturation at a node can only start to increase when the pore pressure at the node below it becomes positive, and thus full saturation is reached there.) This distance can be reduced by increasing the number of zones in the column height. The evolution of nodal pore pressure with time follows a stepwise pattern, more pronounced as the fluid is less compressible. This behavior occurs because a node must be fully saturated before its pore pressure can increase. One way to reduce this effect without changing the grid size is to introduce flow in the unsaturated region (capillary pressure) in the fluid flow formulation.

Figure 1: Location of filling front (\(\widehat{x}\) vs \(\widehat{t}\)) — no gravity (analytical solution = solid line; numerical values = dotted lines).

Figure 2: Location of filling front (\(\widehat{x}\) vs \(\widehat{t}\)) — with gravity (analytical solution = solid line; numerical values = dotted lines).

Data Files

1DFillingPorousRegion-NoGravity.f3dat

model new
model largestrain off
model title "One-dimensional filling - no gravity"
fish automatic-create off
model configure fluid
call 'fishFunctions'
@setup
zone create brick size 1 1 25 point 1 (0.025,0,0) point 2 (0,0.025,0) point 3 (0,0,0.625)
; --- fluid flow model ---
zone fluid cmodel assign isotropic
zone fluid property permeability @c_perm porosity @c_poro biot 0. ; 
zone fluid biot on
zone gridpoint initialize biot 100
zone gridpoint initialize saturation 0.0
zone face apply pore-pressure @c_p0 range position-z 0
; --- settings ---
model mechanical active off
model fluid active on
@solution1
table 1 label 'Analytical solution'
table 2 label ' 1% saturation front'
table 3 label '99% saturation front'
; --- test ---
fish set @lzf = .025 @uzf = .025
model solve fluid time-total 0.25
model save 'asat1'
return

1DFillingPorousRegion-Gravity.f3dat

model new
model largestrain off
model title "One-dimensional filling - with gravity"
fish automatic-create off
model configure fluid
call 'fishFunctions'
@setup
zone create brick size 1 1 25 point 1 (0.025,0,0) point 2 (0,0.025,0) point 3 (0,0,0.625)
; --- fluid flow model ---
zone fluid cmodel assign isotropic
zone fluid property permeability @c_perm poros @c_poro biot 0. 
zone fluid biot on
zone gridpoint initialize biot 100
zone gridpoint initialize saturation 0.0
zone face apply pore-pressure @c_p0 range position-z 0
; --- settings ---
zone initialize fluid-density @c_den
model gravity 0 0 @c_grav
model mechanical active off
model fluid active on
@solution2
table 1 label 'Analytical solution'
table 2 label ' 1% saturation front'
table 3 label '99% saturation front'
; --- test ---
fish set @lzf = .025 @uzf = .025
model solve fluid time-total 0.25
model save 'asat2'
return