This repository compares two approaches to injecting fluid at a fixed flow rate and temperature in MOOSE's PorousFlow Thermo-hydro simulations. The first approach fixes the temperature at the inlet using a DirichletBC, while the second approach injects an amount of enthalpy based on the inlet fluid temperature. This second, enthalpy-based approach can be done in the 'stupid way' using a PorousFlowSink BC or it can be done more effectively with the newly-developed PorousFlowSinkEnthalpy boundary condition. The advantage of PorousFlowSinkEnthalpy is that the inlet temperature can be specified, even when the EOS is complicated, which would not allow an a-priori calcualation of the enthalpy flow rate. The PorousFlowSinkEnthalpy is in the following branch: https://github.com/andrsd/moose/tree/ent-inj-bc. Many thanks to David Andrs for his development of this!
We use the SimpleFluidProperties equation of state: Viscosity = 1 Pa-s, heat capacity = 4200 J/kg/K, reference density = 1000 kg/m3, bulk modulus of the fluid = 2e9 Pa. The rock material properties are: permeability = 10^-7 m2, porosity = 0.1, rock density = 2500 kg/m3, specific heat capacity = 800, dry thermal conductivity = 2.1, wet thermal conductivity = 1.8. The domain is 200 x 1 x 200 m and the area of the injection and extraction BCs are 200 m2.
The domain is initially at zero degrees and has a BCs that encourage a flow of 40000 kg/s, or 40 m3/s, across it. After the initialization period, the inlet temperature is raised in one of three ways: (a) the temperature at the boundary is raised to 1 degree (input_PresetBC.i), (b) the relevant enthalpy injection rate is specified (i.e. the enthalpy from 40000 kg/s at 1 degree, in input_PFSink.i), or (c) in input_enthalpy.i the temperature is specified (using the T_in parameter), which automatically calculates the correct enthalpy of the fluid at the inlet. The outflow boundary produces fluid at 40000 kg/s and allows temperature to leave with the fluid. This is achieved by having very similar PorousFlowSink cards at the outlet for temperature and pressure. The only difference are the variable names (temp instead of pp) and that use_enthalpy is included for the temperature condition.
[./injection]
type = PorousFlowSink
variable = pp
boundary = left
flux_function = -200
fluid_phase = 0
use_mobility = false
save_in = fluxes_in
[../]
[./injection_heat]
type = PresetBC
variable = temp
boundary = left
value = 1.
save_in = heat_fluxes_in
[../]
The enthalpy injection rate is the product of the flow rate, the fluid heat capacity, and the temperature change = 40,000 kg/s * 4200 J/kg/K * 1 K = 1.68e8 J/s. Dividing by the area gives the flux_function = 840000 J/m2/s for the PorousFlowSink BC in input_PFSink.i. Note that this was an a-priori calculation, but it would not be possible for EOS where the enthalpy is a more complex function of pressure or temperature (for example Water97FluidPropertiesFluid, or CO2 or methane).
[./injection]
type = PorousFlowSink
variable = pp
boundary = left
flux_function = -200 ## injeciton of fluid at 200 kg/m2/s
fluid_phase = 0
use_mobility = false
save_in = fluxes_in
[../]
[./injection_heat]
type = PorousFlowSink
variable = temp
boundary = left
flux_function = -840000. #This is the amount of enthalpy calculated a priori
fluid_phase = 0
use_mobility = false
use_enthalpy = false #NOTE this is false since a priori calculation was already done for flux_function
save_in = heat_fluxes_in
[../]
[./injection]
type = PorousFlowSink
variable = pp
boundary = left
flux_function = -200 ## injection of fluid at 200 kg/m2/s
fluid_phase = 0
use_mobility = false
save_in = fluxes_in
[../]
[./injection_heat]
type = PorousFlowSinkEnthalpy
variable = temp
boundary = left
flux_function = -200 #NOTE - flux_function matches the [./injection] card above
pressure = pp
T_in = 1 # Temp of inlet fluid is specified here
fp = simple_fluid
save_in = heat_fluxes_in
[../]
The following BCs are used in all the input files, at the outflow boundary.
[./production]
type = PorousFlowSink
variable = pp
boundary = right
flux_function = 200 ## extraction of fluid at 200 kg/m2/s
fluid_phase = 0
use_mobility = false
save_in = fluxes_out
[../]
#Note that production_heat is identical to production at outflow boundary
#except that variable = temp and use_enthalpy = true
[./production_heat]
type = PorousFlowSink
variable = temp
boundary = right
flux_function = 200
fluid_phase = 0
use_mobility = false
use_enthalpy = true
save_in = heat_fluxes_out
[../]
The figure below shows temperature versus distance along the flow path, x. Results are shown at t = 1, 10, 100, and 350 seconds. The solid lines are the results when an enthalpy injection rate is specified (with PorousFlowSink), the dashed lines are when the inlet enthalpy is calculated from a specified temperature with PorousFlowSinkEnthalpy, and the dotted lines are results when the temperature is directly specified (with PresetBC). The two enthalpy injection approaches match perfectly, which shows that PorousFlowSinkEnthalpy is working correctly. At early times, the enthalpy injection approach has to heat up the element at the injection. By t = 100, the temperature at the inlet is approximately 1 degree C. On the other hand, when the temperature at the inlet is fixed to 1 degree, the temperature gradient is very large at the inlet, essentially providing an unrealistically large amount of energy at early times. By the time the simulation has reached later times (t >=350 s), the differences between the two approaches are minor. It is likely that the enthalpy specification is more accurate than the fixed temperature BC, but the differences may be inconsequential depending on specifics of the problem (mesh size, time scale, length scale, material properties, etc.). PorousFlowSinkEnthalpy appears to be the best approach to specifying BCs of this nature. (Note the graph and its data are in DataAnalysis.xlsx.)
