Point and line sources/sinks

Injecting fluid at specified temperature

When injecting a fluid at specified temperature (also computed as a postprocessor), users can add another Dirac kernel PorousFlowPointEnthalpySourceFromPostprocessor. (Alternately, users can also fix the temperature of an injected fluid using a Dirichlet BC, but this adds/subtracts heat energy to entire nodal volumes of porous material and fluid, so may lead to unacceptable errors from the additional heat energies added/removed.)

Such example is shown here:

[DiracKernels]
  [source]
    type = PorousFlowPointSourceFromPostprocessor
    variable = pressure
    mass_flux = mass_flux_in
    point = '0.5 0.5 0'
  []
  [source_h]
    type = PorousFlowPointEnthalpySourceFromPostprocessor
    variable = temperature
    mass_flux = mass_flux_in
    point = '0.5 0.5 0'
    T_in = T_in
    pressure = pressure
    fp = simple_fluid
  []
[]

[Postprocessors]
  [total_mass]
    type = PorousFlowFluidMass
    execute_on = 'initial timestep_end'
  []

  [total_heat]
    type = PorousFlowHeatEnergy
  []

  [mass_flux_in]
    type = FunctionValuePostprocessor
    function = mass_flux_in_fn
    execute_on = 'initial timestep_end'
  []

  [avg_temp]
    type = ElementAverageValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []

  [T_in]
    type = FunctionValuePostprocessor
    function = T_in_fn
    execute_on = 'initial timestep_end'
  []
[]

The wellbore pressure

The wellbore pressure is an input into Peaceman's formula. For any along the wellbore, the wellbore pressure is defined as Here

• is an input parameter, bottom_p_or_t. It is the pressure at the bottom of the wellbore.

• is the position (a point in 3D) of the bottom of the wellbore. It is defined to be the last point in the point_file.

• is a weight vector pointing downwards (product of fluid density and gravity), unit_weight. This means that will be the pressure at point in the wellbore, due to gravitational head. If these gravitational effects are undesirable, the user may simply specify (unit_weight = '0 0 0').

Although unusual, PorousFlow also allows to represent a temperature, by setting function_of = temperature, which is useful for non-fluid models that contain polyline sources/sinks of heat.

Peaceman's fluid flux

Peaceman writes as (4) Let us discuss each term on the RHS separately. For boreholes that involve heat only (with function_of = temperature) the in the above expression and the discussion below should be replaced by the temperature.

The mobility

Eq. (4) contains , , and , which are the fluid relative permeability, density and viscosity. Hence the term is the mobility, so users should choose to multiply by the mobility (use_mobility = true) when using PorousFlowPeacemanBorehole. Recall that all the multiplicative factors, including mobility, are evaluated at the nodal positions, not the position . This ensures superior numerical convergence.

The character

The in Eq. (4) is called the character of the wellbore. There are two standard choices (note that depends only on the absolute value ).

• With character = 1 then for , and zero otherwise. In this case the wellbore is called a production wellbore, and it acts as a sink, removing fluid from the porespace.

• With character=-1, for , and zero otherwise. In this case the wellbore is called an injection wellbore, and it acts as a source, adding fluid to the porespace.

Generalising, character may be chosen to be different from , which is useful for time-varying borehole strengths. The above two choices are generalised to:

• if the wellbore is active only for (and has zero strength otherwise);

• if the wellbore is active only for .

The well constant

The function is called the well_constant of the wellbore, and is measured in units of length. Usually it is computed automatically by PorousFlow using Peaceman's formulae, below.

Peaceman described the form of by performing 2D analytical and numerical calculations to obtain the following formula

In this formula:

• the borehole is oriented along the direction (this is generalised in PorousFlow: the direction is defined by the plaintext file of points mentioned in Eq. (2));

• and are the diagonal components of the permeability tensor in the plane (defined by a PorousFlow Material, and generalised in PorousFlow to consider the components normal to the borehole direction);

• is the length of the borehole (defined by the plaintext file of points in Eq. (2));

• is the borehole radius (an input parameter, which is encoded in the plaintext file of points as the first entry in each row: see Eq. (2));

• is the effective borehole radius.

For a cell-centred finite-difference approach, Peaceman (1983) found that Here and are the finite-difference spatial sizes.

Other authors have generalised Peaceman's approach to writing for different geometrical situations. Some of these are contained in Chen and Zhang (2009), where they show that for a finite element situation with square elements of size , the borehole at a nodal position, and isotropic permeability

Note that the 0.113 is substantially different to Peaceman's 0.2, demonstrating that this method of introducing wellbores dependent on the geometry. The user may specify this quantity via the re_constant input parameter (which defaults to Peaceman's 0.28).

Simple examples

In the following examples, a vertical borehole is placed through the centre of a single element, and fluid flow to the borehole as a function of porepressure is measured. The borehole geometry is

0.1 0 0 0.5
0.1 0 0 -0.5

meaning that the borehole radius is 0.1m and it is 1m long.

These examples are part of the MOOSE test suite, and they are testing that Eq. (4) is encoded correctly. The parameters common to each of these tests are:

It is remotely possible that the MOOSE implementation applies the borehole flux incorrectly, but records it as a Postprocessor correctly as specified by Eq. (4). Therefore, these simulations also record the fluid mass and mass-balance error in order to check that the fluid mass is indeed being correctly changed by the borehole.

A production wellbore:

[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole

    # Because the Variable for this Sink is pp, and pp is associated
    # with the fluid-mass conservation equation, this sink is extracting
    # fluid mass (and not heat energy or something else)
    variable = pp

    # The following specfies that the total fluid mass coming out of
    # the porespace via this sink in this timestep should be recorded
    # in the pls_total_outflow_mass UserObject
    SumQuantityUO = borehole_total_outflow_mass

    # The following file defines the polyline geometry
    # which is just two points in this particular example
    point_file = bh02.bh

    # First, we want Peacemans f to be a function of porepressure (and not
    # temperature or something else).  So bottom_p_or_t is actually porepressure
    function_of = pressure
    fluid_phase = 0

    # The bottomhole pressure
    bottom_p_or_t = 0

    # In this example there is no increase of the wellbore pressure
    # due to gravity:
    unit_weight = '0 0 0'

    # PeacemanBoreholes should almost always have use_mobility = true
    use_mobility = true

    # This is a production wellbore (a sink of fluid that removes fluid from porespace)
    character = 1
  []
[]

with an initially fully-saturated medium yields the correct result (Figure 1 and Figure 2):

Figure 1: The flow to a production wellbore in MOOSE agrees with the expected result.

Figure 2: The mass-balance error is virtually zero.

An injection wellbore:

[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = pp
    SumQuantityUO = borehole_total_outflow_mass
    point_file = bh03.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 'insitu_pp'
    unit_weight = '0 0 0'
    use_mobility = true
    character = -1
  []
[]

with a fully-saturated medium yields the correct result (Figure 3 and Figure 4):

Figure 3: The flow from an injection wellbore in MOOSE agrees with the expected result.

Figure 4: The mass-balance error is virtually zero.

A production wellbore with bottomhole porepressure MPa so it forces desaturation of an initially-saturated medium:

[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = pp
    SumQuantityUO = borehole_total_outflow_mass
    point_file = bh02.bh
    fluid_phase = 0
    bottom_p_or_t = -1E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
[]

yields the correct result (Figure 5 and Figure 6):

Figure 5: The flow to a production wellbore in MOOSE agrees with the expected result.

Figure 6: The mass-balance error is virtually zero.

An injection wellbore with so it is only active when the rock porepressure is negative:

[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = pp
    SumQuantityUO = borehole_total_outflow_mass
    point_file = bh03.bh
    fluid_phase = 0
    bottom_p_or_t = 0
    unit_weight = '0 0 0'
    use_mobility = true
    character = -1
  []
[]

yields the correct result for an initially unsaturated medium (Figure 7 and Figure 8):

Figure 7: The flow from an injection wellbore in MOOSE agrees with the expected result.

Figure 8: The mass-balance error is virtually zero.

Each of these record the total fluid flux (kg) injected by or produced by the borehole in a PorousFlowSumQuantity UserObject and outputs this result using a PorousFlowPlotQuantity Postprocessor:

  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []

  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []

Reproducing the steady-state 2D analytical solution

The PorousFlow fluid equation (see governing equations) for a fully-saturated medium with and large constant bulk modulus becomes Darcy's equation where , with notation described in the nomenclature. In the isotropic case (where ), the steadystate equation is just Laplace's equation

Place a borehole of radius and infinite length oriented along the axis. Then the situation becomes 2D and can be solved in cylindrical coordinates, with and independent of . If the pressure at the borehole wall is , then the fluid density is . Assume that at the fluid pressure is held fixed at , or equivalently the density is held fixed at . Then the solution of Laplace's equation is well-known to be (5) This is the fundamental solution used by Peaceman and others to derive expressions for by comparing with numerical expressions resulting from Eq. (4).

Chen and Zhang Chen and Zhang (2009) have derived an expression for in the case where this borehole is placed at a node in a square mesh. This test compares the MOOSE steadystate solution with a single borehole with defined by Chen and Zhang's formula (specifically, the re_constant needs to be changed from its default value) is compared with Eq. (5) to illustrate that the MOOSE implementation of a borehole is correct.

The PorousFlowPeacemanBorehole is:

[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = pp
    SumQuantityUO = borehole_total_outflow_mass
    point_file = bh07.bh
    fluid_phase = 0
    bottom_p_or_t = 0
    unit_weight = '0 0 0'
    use_mobility = true
    re_constant = 0.1594 # use Chen and Zhang version
    character = 2 # double the strength because bh07.bh only fills half the mesh
  []
[]

The mesh used is shown in Figure 9.

Figure 9: The mesh used in the comparison with Eq. (5), with the green dot indicating the position of the borehole. The central elements are m, and the outer boundary is at radius 300m.

Figure 10 show the comparison between the MOOSE result and Eq. (5). Most parameters in this study are identical to those given in the Table 1 with the following exceptions: the mesh is shown in Figure 9; the permeability is m; the borehole radius is 1 m; the borehole pressure is ; the outer radius is m; and the outer pressure is MPa.

Figure 10: Comparison of the MOOSE results (dots) with the analytical solution Eq. (5) for the steadystate porepressure distribution surrounding single borehole.

Injecting and producing in thermo-hydro simulations

The PorousFlowPeacemanBorehole may be used to model injection and production in thermo-hydro simulations. Suppose that cool water is being injected into an initially hot reservoir via a vertical injection well, and that water is being removed from the system via a vertical production well. In this example we shall study an essentially 2D situation without gravity with parameters defined in Table 2

The fluid properties are defined as:

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    thermal_expansion = 2E-4
    bulk_modulus = 2E9
    viscosity = 1E-3
    density0 = 1000
    cv = 4000.0
    cp = 4000.0
  []
[]

The fluid injection and production is implemented in a way that is now familiar:

  [fluid_injection]
    type = PorousFlowPeacemanBorehole
    variable = porepressure
    SumQuantityUO = injected_mass
    point_file = injection.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 21E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = -1
  []
  [fluid_production]
    type = PorousFlowPeacemanBorehole
    variable = porepressure
    SumQuantityUO = produced_mass
    point_file = production.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []

The injection temperature is set via:

[BCs]
  [injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 300
    boundary = central_nodes
  []
[]

Alternatively, the heat energy of the injected fluid could be worked out and injected into the heat-conservation equation via a borehole.

The production of fluid means that heat energy must be removed at exactly the rate associated with the fluid-mass removal. This is implemented by using an identical PorousFlowPeacemanBorehole to the fluid production situation but associating it with the temperature variable and with use_enthalpy = true:

  [remove_heat_at_production_well]
    type = PorousFlowPeacemanBorehole
    variable = temperature
    SumQuantityUO = produced_heat
    point_file = production.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    use_enthalpy = true
    character = 1

Here the heat energy per timestep is saved into a PorousFlowSumQuantity UserObject which may be extracted using a PorousFlowPlotQuantity Postprocessor, which would be an important quantity in a geothermal application:

[Postprocessors]
  [heat_joules_extracted_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_heat
  []
[]

Scaling of the variables ensures good convergence:

[Variables]
  [porepressure]
    initial_condition = 20E6
  []
  [temperature]
    initial_condition = 400
    scaling = 1E-6 # fluid enthalpy is roughly 1E6
  []
[]

Some results (run on a fine mesh) are shown in Figure 11 and Figure 12.

Figure 11: Porepressure after s of injection and production. The black spots indicate the wells.

Figure 12: Temperature after s of injection and production. the black spots indicate the wells.