Injection ala Buckley and Leverett

PorousFlow is compared with a simple single-phase one-dimensional Buckley-Leverett problem (Buckley and Leverett, 1942). The single-phase fluid moves in a region m. A fully-saturated front initially sits at position , while the region is initially unsaturated. With zero suction function , there is no diffusion of the sharp front, and it progresses towards the right (increasing ).

# Buckley-Leverett 1-phase.
# The front starts at (around) x=5, and at t=50 it should
# have moved to x=9.6.  The version below has a nonzero
# suction function, and at t=50, the front sits between
# (about) x=9.6 and x=9.9.  Changing the van-Genuchten
# al parameter to 1E-4 softens the front so it sits between
# (about) x=9.7 and x=10.4, and the simulation runs much faster.
# With al=1E-2 and nx=600, the front sits between x=9.6 and x=9.8,
# but takes about 100 times longer to run.

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 1
  nx = 150
  xmin = 0
  xmax = 15
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
  compute_enthalpy = false
  compute_internal_energy = false
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
    [InitialCondition]
      type = FunctionIC
      function = 'max((1000000-x/5*1000000)-20000,-20000)'
    []
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
  [flux0]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [left]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 980000
  []
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [sat]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../../../../syntax/AuxKernels/index.html"}>>>]
  [sat]
    type = MaterialStdVectorAux<<<{"description": "Extracts a component of a material type std::vector<Real> to an aux variable.  If the std::vector is not of sufficient size then zero is returned", "href": "../../../../source/auxkernels/MaterialStdVectorAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = sat
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
    index<<<{"description": "The index to consider for this kernel"}>>> = 0
    property<<<{"description": "The material property name."}>>> = PorousFlow_saturation_qp
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.8
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-3
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e6
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pp
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-10 0 0  0 1E-10 0  0 0 1E-10'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.15
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  active<<<{"description": "If specified only the blocks named will be visited and made active"}>>> = andy
  [andy]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'gmres bjacobi 1E-10 1E-10 20'
  []
[]

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [timestepper]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../../../../source/functions/PiecewiseLinear.html"}>>>
    x<<<{"description": "The abscissa values"}>>> = '0    0.01 0.1 1   1.5 2   20  30  40  50'
    y<<<{"description": "The ordinate values"}>>> = '0.01 0.1  0.2 0.3 0.1 0.3 0.3 0.4 0.4 0.5'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  end_time = 50
  [TimeStepper<<<{"href": "../../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = FunctionDT
    function = timestepper
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [pp]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0 0'
    end_point<<<{"description": "The ending of the line"}>>> = '15 0 0'
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 150
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = pp
  []
  [sat]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    warn_discontinuous_face_values<<<{"description": "Whether to return a warning if a discontinuous variable is sampled on a face"}>>> = false
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0 0'
    end_point<<<{"description": "The ending of the line"}>>> = '15 0 0'
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 150
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = sat
  []
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = bl01
  [csv]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    sync_only<<<{"description": "Only export results at sync times"}>>> = true
    sync_times<<<{"description": "Times at which the output and solution is forced to occur"}>>> = '0.01 50'
  []
  [exodus]
    type = Exodus<<<{"description": "Object for output data in the Exodus format", "href": "../../../../source/outputs/Exodus.html"}>>>
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'initial final'
  []
[]
(modules/porous_flow/test/tests/buckley_leverett/bl01.i)

This is a difficult problem to simulate numerically as maintaining the sharp front is hard. The front's speed is independent of the relative permeability, since the fluid is flowing from a fully-saturated region (where ). This problem is therefore a good test of the full upwinding.

In the simulation listed above, the pressure at the left boundary is kept fixed at MPa, while the right boundary is kept fixed at Pa, so the difference is MPa. The medium's permeability is set to and its porosity is . It is not possible to use a zero suction function in the MOOSE implementation, but using the van Genuchten parameters Pa and approximates it. The fluid viscosity is Pa.s.

The initial condition is which is shown in Figure 1. With the suction function defined above this gives

Figure 1: Initial setup of the Buckley-Leverett problem where porepressure is a piecewise linear function. The region on the left is fully saturated, while the region on the right has saturation 0.000006. During simulation the value of porepressure on the left boundary is held fixed.

Good approximations for the pressure and the front position may be determined from (1) which has solution For the parameters listed above, the front will be at position m at s. This solution is only valid for zero capillary suction. A nonzero suction function will tend to diffuse the sharp front.

With coarse meshes it is impossible to simulate a sharp front, of course, since the front is spread over at least one element.

Figure 2 shows the results from a MOOSE simulation with a uniform mesh of size m. At s the front in this simulation sits at m as desired. However, the simulation takes seconds to complete due to the very low values of saturation obtained for van Genuchten Pa. Other simulations give similar results but run much faster. For instance, the test suite listed above uses the van Genuchten parameter Pa, but the front diffuses a little into the unsaturated region (the front sits between m and m at s).

Figure 2: The MOOSE solution of the Buckley-Leverett problem at time 50s. Left: saturation. Right: porepressure. The front sits at x 9.6m as expected from the analytical solution.

References

  1. E. Buckley and M. Leverett. Mechanism of fluid displacement in sands. Transactions of the AIME, 146:107–116, 04 1942. doi:10.2118/942107-G.[BibTeX]