# Heat energy produced by plastic deformation

In PorousFlow, it is assumed that plastic deformation causes a heat energy-density rate (J.m.s) of (1) where is a coefficient (s), is the porosity, is the stress, and is the plastic strain. This is implemented in the PorousFlowPlasticHeatEnergy Kernel

In the set of simple tests described below, the porous skeleton contains no fluid, and no heat flow is considered, so the heat energy released simply heats up the porous skeleton: (2) The porosity () and volumetric heat capacity of the rock grains () are chosen to be constant.

Perfect, capped, weak-plane plasticity is used, so that the admissible zone is defined by (3) Here is the tensile strength, is the compressive strength, is the cohesion and is the friction angle. The elastic tensor is chosen to be isotropic (4) The parameters in these expressions are chosen to be: , , , , , and (all in consistent units).

In each experiment a single finite-element of unit size is used.

## Tensile failure

The top of the finite element is pulled upwards with displacement: (5) while the displacement on the element's bottom, and in the and directions is chosen to be zero.

# Tensile heating, using capped weak-plane plasticity
# z_disp(z=1) = t
# totalstrain_zz = t
# with C_ijkl = 0.5 0.25
# stress_zz = t, but with tensile_strength = 1, stress_zz = min(t, 1)
# so plasticstrain_zz = t - 1
# heat_energy_rate = coeff * (t - 1)
# Heat capacity of rock = specific_heat_cap * density = 4
# So temperature of rock should be:
# (1 - porosity) * 4 * T = (1 - porosity) * coeff * (t - 1)
[Mesh]
type = GeneratedMesh
dim = 3
xmin = -10
xmax = 10
zmin = 0
zmax = 1
[]

[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
[]

[Variables]
[./temperature]
[../]
[]

[Kernels]
[./energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = temperature
[../]
[./phe]
type = PorousFlowPlasticHeatEnergy
variable = temperature
[../]
[]

[AuxVariables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]

[AuxKernels]
[./disp_z]
type = FunctionAux
variable = disp_z
function = z*t
[../]
[]

[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = temperature
number_fluid_phases = 0
number_fluid_components = 0
[../]
[./coh]
type = TensorMechanicsHardeningConstant
value = 100
[../]
[./tanphi]
type = TensorMechanicsHardeningConstant
value = 1.0
[../]
[./t_strength]
type = TensorMechanicsHardeningConstant
value = 1
[../]
[./c_strength]
type = TensorMechanicsHardeningConstant
value = 1
[../]
[]

[Materials]
[./rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2
density = 2
[../]
[./temp]
type = PorousFlowTemperature
temperature = temperature
[../]
[./porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[../]
[./phe]
type = ComputePlasticHeatEnergy
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0.5 0.25'
[../]
[./strain]
type = ComputeIncrementalSmallStrain
displacements = 'disp_x disp_y disp_z'
[../]
type = ComputeMultipleInelasticStress
inelastic_models = mc
perform_finite_strain_rotations = false
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanphi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-10
perfect_guess = true
[../]
[]

[Postprocessors]
[./temp]
type = PointValue
point = '0 0 0'
variable = temperature
[../]
[]

[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000'
[../]
[]

[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 10
[]

[Outputs]
file_base = tensile01
csv = true
[]

(modules/porous_flow/test/tests/plastic_heating/tensile01.i)

This implies the only non-zero component of total strain is (6) The constitutive law implies (7) This stress is admissible for , while for the system yields in tension: (8) and the plastic strain is, (9) This means that the material's temperature should increase as (10) with no temperature incraese for for .

This result is reproduced exactly by PorousFlow.

## Compressive failure

The top of the finite element is pushed downwards with displacement: (11) while the displacement on the element's bottom, and in the and directions is chosen to be zero.

# Tensile heating, using capped weak-plane plasticity
# z_disp(z=1) = -t
# totalstrain_zz = -t
# with C_ijkl = 0.5 0.25
# stress_zz = -t, but with compressive_strength = 1, stress_zz = max(-t, -1)
# so plasticstrain_zz = -(t - 1)
# heat_energy_rate = coeff * (t - 1)
# Heat capacity of rock = specific_heat_cap * density = 4
# So temperature of rock should be:
# (1 - porosity) * 4 * T = (1 - porosity) * coeff * (t - 1)
[Mesh]
type = GeneratedMesh
dim = 3
xmin = -10
xmax = 10
zmin = 0
zmax = 1
[]

[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
[]

[Variables]
[./temperature]
[../]
[]

[Kernels]
[./energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = temperature
[../]
[./phe]
type = PorousFlowPlasticHeatEnergy
variable = temperature
[../]
[]

[AuxVariables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]

[AuxKernels]
[./disp_z]
type = FunctionAux
variable = disp_z
function = '-z*t'
[../]
[]

[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = temperature
number_fluid_phases = 0
number_fluid_components = 0
[../]
[./coh]
type = TensorMechanicsHardeningConstant
value = 100
[../]
[./tanphi]
type = TensorMechanicsHardeningConstant
value = 1.0
[../]
[./t_strength]
type = TensorMechanicsHardeningConstant
value = 1
[../]
[./c_strength]
type = TensorMechanicsHardeningConstant
value = 1
[../]
[]

[Materials]
[./rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2
density = 2
[../]
[./temp]
type = PorousFlowTemperature
temperature = temperature
[../]
[./porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[../]
[./phe]
type = ComputePlasticHeatEnergy
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0.5 0.25'
[../]
[./strain]
type = ComputeIncrementalSmallStrain
displacements = 'disp_x disp_y disp_z'
[../]
type = ComputeMultipleInelasticStress
inelastic_models = mc
perform_finite_strain_rotations = false
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanphi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-10
perfect_guess = true
[../]
[]

[Postprocessors]
[./temp]
type = PointValue
point = '0 0 0'
variable = temperature
[../]
[]

[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000'
[../]
[]

[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 10
[]

[Outputs]
file_base = compressive01
csv = true
[]

(modules/porous_flow/test/tests/plastic_heating/compressive01.i)

This implies only non-zero component of total strain is (12) The constitutive law implies (13) This stress is admissible for , while for the system yields in compression (14) and the plastic strain is, (15) This means that the material's temperature should increase as (16) and there should be no temperature increase for .

MOOSE produces this result exactly

## Shear failure

The top of the finite element is sheared with displacement: (17) while the displacement on the element's bottom, and in the and directions is chosen to be zero.

# Tensile heating, using capped weak-plane plasticity
# x_disp(z=1) = t
# totalstrain_xz = t
# with C_ijkl = 0.5 0.25
# stress_zx = stress_xz = 0.25*t, so q=0.25*t, but
# with cohesion=1 and tan(phi)=1: max(q)=1.  With tan(psi)=0,
# the plastic return is always to (p, q) = (0, 1),
# so plasticstrain_zx = max(t - 4, 0)
# heat_energy_rate = coeff * (t - 4) for t>4
# Heat capacity of rock = specific_heat_cap * density = 4
# So temperature of rock should be:
# (1 - porosity) * 4 * T = (1 - porosity) * coeff * (t - 4)
[Mesh]
type = GeneratedMesh
dim = 3
xmin = -10
xmax = 10
zmin = 0
zmax = 1
[]

[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
[]

[Variables]
[./temperature]
[../]
[]

[Kernels]
[./energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = temperature
[../]
[./phe]
type = PorousFlowPlasticHeatEnergy
variable = temperature
coeff = 8
[../]
[]

[AuxVariables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]

[AuxKernels]
[./disp_x]
type = FunctionAux
variable = disp_x
function = 'z*t'
[../]
[]

[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = temperature
number_fluid_phases = 0
number_fluid_components = 0
[../]
[./coh]
type = TensorMechanicsHardeningConstant
value = 1
[../]
[./tanphi]
type = TensorMechanicsHardeningConstant
value = 1.0
[../]
[./tanpsi]
type = TensorMechanicsHardeningConstant
value = 0.0
[../]
[./t_strength]
type = TensorMechanicsHardeningConstant
value = 10
[../]
[./c_strength]
type = TensorMechanicsHardeningConstant
value = 10
[../]
[]

[Materials]
[./rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2
density = 2
[../]
[./temp]
type = PorousFlowTemperature
temperature = temperature
[../]
[./porosity]
type = PorousFlowPorosityConst
porosity = 0.7
[../]
[./phe]
type = ComputePlasticHeatEnergy
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0.5 0.25'
[../]
[./strain]
type = ComputeIncrementalSmallStrain
displacements = 'disp_x disp_y disp_z'
[../]
type = ComputeMultipleInelasticStress
inelastic_models = mc
perform_finite_strain_rotations = false
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-10
perfect_guess = true
[../]
[]

[Postprocessors]
[./temp]
type = PointValue
point = '0 0 0'
variable = temperature
[../]
[]

[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000'
[../]
[]

[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 10
[]

[Outputs]
file_base = shear01
csv = true
[]

(modules/porous_flow/test/tests/plastic_heating/shear01.i)

This implies only non-zero component of total strain is (18) The constitutive law implies (19) This stress is admissible for , while for the system yields in shear: (20) and the plastic strain is, (21) This means that the material's temperature should increase as (22) while the right-hand side is zero for .

PorousFlow produces this result exactly.