- PorousFlowDictatorThe UserObject that holds the list of PorousFlow variable names
C++ Type:UserObjectName
Controllable:No
Description:The UserObject that holds the list of PorousFlow variable names
- matrix_reference_temperatureUndisturbed temperature of the matrix adjacent to the fracture. This should usually be a linear-lagrange variable or a constant
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Undisturbed temperature of the matrix adjacent to the fracture. This should usually be a linear-lagrange variable or a constant
- matrix_thermal_diffusivityThermal diffusivity of the matrix adjacent to the fracture (thermal conductivity divided by volumetric heat capacity), m^2/s. This is the diffusivity in the direction normal to the fracture (or the matrix is assumed isotropic), and it is assumed constant
C++ Type:Real
Unit:(no unit assumed)
Range:matrix_thermal_diffusivity > 0
Controllable:No
Description:Thermal diffusivity of the matrix adjacent to the fracture (thermal conductivity divided by volumetric heat capacity), m^2/s. This is the diffusivity in the direction normal to the fracture (or the matrix is assumed isotropic), and it is assumed constant
- porosity_refThe porosity at reference effective porepressure, temperature and volumetric strain. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The porosity at reference effective porepressure, temperature and volumetric strain. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)
- thermal_opening_coeffChange in porosity per unit heat content of the adjacent matrix. The heat content, I = int (T - matrix_reference_temperature) dz, has units K.m, so when porosity represents the aperture of a lower-dimensional fracture this coefficient has units 1/K. Porosity changes by -thermal_opening_coeff * I, so cooling (I < 0) opens the fracture when this coefficient is positive
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Change in porosity per unit heat content of the adjacent matrix. The heat content, I = int (T - matrix_reference_temperature) dz, has units K.m, so when porosity represents the aperture of a lower-dimensional fracture this coefficient has units 1/K. Porosity changes by -thermal_opening_coeff * I, so cooling (I < 0) opens the fracture when this coefficient is positive
PorousFlowPorosityVinsomeWesterveld
This Material was contributed by Bernard Che Ngu (Baylor University).
When cold fluid is injected into a fractured rock, the rock next to each fracture cools and contracts, and the fracture opens. This Material models that effect. It is intended for models in which a fracture is represented by lower-dimensional elements, either in the same mesh as the matrix or in a separate MultiApp (see flow through fractured media), and the porosity of those elements is the fracture aperture (multiplied by the fracture porosity, which is usually unity). The Material accounts for conductive heat transfer between the fracture and the matrix. It does not account for fluid flow between the fracture and the matrix, or the heat carried by such flow.
The Material computes porosity (at the nodes or quadrature points, depending on the "at_nodes" flag) according to
(1)
where is the porosity given by PorousFlowPorosityLinear (see that page for the meaning of the other symbols), is the "thermal_opening_coeff", and
(2)
is the heat content (in Km) of the matrix adjacent to the fracture. Here is the distance in the normal direction from the fracture face into the matrix and is the "matrix_reference_temperature", which is the undisturbed matrix temperature. When the matrix cools, , so a positive increases the porosity (opens the fracture).
The thermal opening depends on how much of the matrix has cooled, not only on the temperature at the fracture face. The cooled layer is often only centimetres to metres thick, which is much thinner than typical matrix elements, so it is not resolved by the mesh (see also the discussion of "retardation" by the matrix in flow through fractured media). Therefore is computed using a semi-analytic temperature profile, as described below.
The Vinsome-Westerveld method
Following Vinsome and Westerveld (1980), who developed the method to compute heat losses to cap and base rock in thermal reservoir simulators (including backflow of heat and cyclic temperature changes), the matrix temperature profile normal to the fracture is approximated by
(3)
where is the temperature change at the fracture face, is the "matrix_thermal_diffusivity" (thermal conductivity divided by volumetric heat capacity) in the direction normal to the fracture, and is the time since "vw_start_time". The exponential factor captures the diffusive penetration of the temperature change into the matrix, and and adjust the shape of the profile as the history of evolves. The quantities , , and depend on but not on , so all of the -dependence of the profile is shown explicitly in Eq. (3). With this profile, the heat content is
The coefficients and are found at each time step by requiring that (1) the heat-conduction equation holds at the fracture face, and (2) energy is conserved, that is, the change in over the time step equals the conductive heat flux through the fracture face, . Discretising these with an implicit (backward-Euler) time step of size and eliminating gives
(4)
where the superscript "old" denotes values at the previous time step. and are stored as stateful material properties, and are both zero initially, which means the matrix is at at the start of the simulation.
Because and are fixed during a time step, Eq. (4) shows that is a linear function of the current temperature. The derivative is therefore a single number at each node, and the Jacobian contribution of the thermal-opening term is as simple to compute as that of the temperature term in PorousFlowPorosityLinear. Vinsome and Westerveld (1980) note that this linearity makes the method easy to couple implicitly to the temperature, which they found to be stable at large time steps, whereas explicit coupling was not.
Choosing the thermal-opening coefficient
If the matrix on both sides of the fracture contracts freely with linear thermal expansion coefficient , each fracture face moves by , so
If the matrix is instead constrained laterally (zero strain parallel to the fracture, and zero stress normal to it), the normal thermal strain is increased by the factor , where is the drained Poisson's ratio, so
These values assume the fracture opening is not resisted by the in-situ stress, so they are upper bounds. In practice, is often treated as a calibration parameter. If the porosity of the lower-dimensional elements is the aperture multiplied by a fracture porosity , then should also be multiplied by .
Assumptions and limitations
Heat conduction in the matrix is one-dimensional, normal to the fracture, and the matrix is semi-infinite. The thermal diffusivity is the value in the direction normal to the fracture (or the matrix is assumed isotropic), and it is constant: it does not depend on temperature or saturation.
Only conductive heat transfer between the fracture and the matrix is considered. Fluid flow between the fracture and the matrix is not included in the calculation of . When the cooled zones of neighbouring fractures overlap, will be overestimated in magnitude.
The matrix is at "matrix_reference_temperature" at "vw_start_time". Before "vw_start_time" the thermal-opening term is zero, and any difference between the fracture-face temperature and "matrix_reference_temperature" at "vw_start_time" is treated as a step change at that time. Because the diffusion length is measured from "vw_start_time", this should be the time at which the fracture-face temperature first starts to change (for instance, when injection begins), as in Vinsome and Westerveld (1980).
The profile Eq. (3) is an approximation. In the example below, for a step change in the fracture-face temperature, the computed is smaller in magnitude than the exact value once the time step is small compared with the elapsed time, and the error is larger when the time step is comparable to the elapsed time (see below).
The heat lost by the matrix is not removed from the simulation by this Material: the thermal-opening term only modifies porosity. The heat exchange between the fracture and the matrix is computed by the usual PorousFlow kernels.
Verification
For a step change in the fracture-face temperature at , the exact matrix temperature is , and integrating over gives . The Vinsome-Westerveld equations can also be solved exactly in this case, if the time step is infinitesimally small. Writing and , the heat-conduction equation at (where because is constant) gives , so . The energy balance , with , then gives and . The ratio of the two is , so the assumed profile underestimates the magnitude of by at all times.
With a finite time step there is an additional error. For a step change, this depends only on the number of time steps since the step (), not on or on the size of , as shown in Table 1.
Table 1: Error in the matrix heat content for a step change in the fracture-face temperature
| Number of time steps since the step | 1 | 10 | 20 | 40 | 100 | 1000 | 10000 |
|---|---|---|---|---|---|---|---|
| Error in |
The input file for the step-change test (K, ms) is
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[porosity_VW]
type = PorousFlowPorosityVinsomeWesterveld<<<{"description": "This Material calculates porosity (usually the aperture of a lower-dimensional fracture) using the PorousFlowPorosityLinear relationship plus a thermal-opening term proportional to the heat content of the adjacent matrix, which is computed using the semi-analytic method of Vinsome and Westerveld (1980)", "href": "PorousFlowPorosityVinsomeWesterveld.html"}>>>
porosity_ref<<<{"description": "The porosity at reference effective porepressure, temperature and volumetric strain. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)"}>>> = 0.1
thermal_opening_coeff<<<{"description": "Change in porosity per unit heat content of the adjacent matrix. The heat content, I = int (T - matrix_reference_temperature) dz, has units K.m, so when porosity represents the aperture of a lower-dimensional fracture this coefficient has units 1/K. Porosity changes by -thermal_opening_coeff * I, so cooling (I < 0) opens the fracture when this coefficient is positive"}>>> = 1E-5
matrix_thermal_diffusivity<<<{"description": "Thermal diffusivity of the matrix adjacent to the fracture (thermal conductivity divided by volumetric heat capacity), m^2/s. This is the diffusivity in the direction normal to the fracture (or the matrix is assumed isotropic), and it is assumed constant"}>>> = 1E-6
matrix_reference_temperature<<<{"description": "Undisturbed temperature of the matrix adjacent to the fracture. This should usually be a linear-lagrange variable or a constant"}>>> = 200
[]
[]Descriptions of other porosity classes can be found in the porosity documentation
Input Parameters
- P_coeff0Effective porepressure coefficient
Default:0
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Effective porepressure coefficient
- P_refThe reference effective porepressure. This should usually be a linear-lagrange variable
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The reference effective porepressure. This should usually be a linear-lagrange variable
- T_coeff0Temperature coefficient
Default:0
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Temperature coefficient
- T_refThe reference temperature. This should usually be a linear-lagrange variable
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The reference temperature. This should usually be a linear-lagrange variable
- at_nodesFalseEvaluate Material properties at nodes instead of quadpoints
Default:False
C++ Type:bool
Controllable:No
Description:Evaluate Material properties at nodes instead of quadpoints
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
- constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
- declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
- epv_coeff0Volumetric-strain coefficient
Default:0
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Volumetric-strain coefficient
- epv_refThe reference volumetric strain. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The reference volumetric strain. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)
- porosity_min0Minimum allowed value of the porosity: if the computed porosity is less than this value, porosity is set to this value instead. The default of zero means porosity is never negative. A negative value may be used, for instance, when comparing with analytical solutions of linear poroelasticity in which porosity can become negative
Default:0
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Minimum allowed value of the porosity: if the computed porosity is less than this value, porosity is set to this value instead. The default of zero means porosity is never negative. A negative value may be used, for instance, when comparing with analytical solutions of linear poroelasticity in which porosity can become negative
- strain_at_nearest_qpFalseWhen calculating nodal porosity that depends on strain, use the strain at the nearest quadpoint. This adds a small extra computational burden, and is not necessary for simulations involving only linear lagrange elements. If you set this to true, you will also want to set the same parameter to true for related Kernels and Materials
Default:False
C++ Type:bool
Controllable:No
Description:When calculating nodal porosity that depends on strain, use the strain at the nearest quadpoint. This adds a small extra computational burden, and is not necessary for simulations involving only linear lagrange elements. If you set this to true, you will also want to set the same parameter to true for related Kernels and Materials
- vw_start_time0Time at which the temperature at the fracture face first starts to change (for instance, when injection begins). The matrix adjacent to the fracture is at matrix_reference_temperature at this time, and the thermal-opening term is zero before it
Default:0
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Time at which the temperature at the fracture face first starts to change (for instance, when injection begins). The matrix adjacent to the fracture is at matrix_reference_temperature at this time, and the thermal-opening term is zero before it
Optional Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- search_methodnearest_node_connected_sidesChoice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
Default:nearest_node_connected_sides
C++ Type:MooseEnum
Controllable:No
Description:Choice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
- zero_modifier0.001If the porosity_min floor is active, the porosity derivatives are set to zero_modifier times their unfloored values (rather than exactly zero) to hint to the nonlinear solver that porosity is not strictly constant, which aids convergence
Default:0.001
C++ Type:Real
Unit:(no unit assumed)
Range:zero_modifier >= 0
Controllable:No
Description:If the porosity_min floor is active, the porosity derivatives are set to zero_modifier times their unfloored values (rather than exactly zero) to hint to the nonlinear solver that porosity is not strictly constant, which aids convergence
Advanced Parameters
- output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
- outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
Outputs Parameters
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Material Property Retrieval Parameters
Input Files
References
- P. K. W. Vinsome and J. Westerveld.
A simple method for predicting cap and base rock heat losses in thermal reservoir simulators.
Journal of Canadian Petroleum Technology, 19(3):87–90, 1980.
doi:10.2118/80-03-04.[Export]