Heat energy produced by plastic deformation

In PorousFlow, it is assumed that plastic deformation causes a heat energy-density rate (J.m $var element = document.getElementById("moose-equation-9dbaa43a-3afd-4c4b-a193-6fd28c2a0841");katex.render("^{-3}", element, {displayMode:false,throwOnError:false});$ .s $var element = document.getElementById("moose-equation-032c67d9-0e47-466f-98e6-820d1fb2b7e4");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ) of $var element = document.getElementById("moose-equation-2e2a78bb-0f82-4cd8-b96c-7dd11a4f2a82");katex.render("c (1-\\phi) \\sigma_{ij}\\epsilon^{\\mathrm{plastic}}_{ij} \\ ,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-b8767d73-34b5-4ca9-b5ad-14ec2fe7e7d9");katex.render("c", element, {displayMode:false,throwOnError:false});$ is a coefficient (s $var element = document.getElementById("moose-equation-ba9b07db-b1ef-499d-be2c-520c7121a5d7");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ), $var element = document.getElementById("moose-equation-47c34973-de26-46ee-9efd-bcf29fd6766f");katex.render("\\phi", element, {displayMode:false,throwOnError:false});$ is the porosity, $var element = document.getElementById("moose-equation-fa42c377-01dc-4339-abb8-d8d2263f5a44");katex.render("\\sigma", element, {displayMode:false,throwOnError:false});$ is the stress, and $var element = document.getElementById("moose-equation-de8f1a10-5dac-4925-8a85-65d4e06a163c");katex.render("\\epsilon^{\\mathrm{plastic}}", element, {displayMode:false,throwOnError:false});$ 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: $var element = document.getElementById("moose-equation-0028ec19-a61f-4387-9974-718de3000d43");katex.render("\\frac{\\partial}{\\partial t} (1 - \\phi)c_{R}\\rho_{R} T = c (1-\\phi) \\sigma_{ij}\\epsilon^{\\mathrm{plastic}}_{ij} \\ .", element, {displayMode:true,throwOnError:false});$ The porosity ( $var element = document.getElementById("moose-equation-e1fb8590-e5b9-4715-8a59-b5a2340c9787");katex.render("\\phi", element, {displayMode:false,throwOnError:false});$ ) and volumetric heat capacity of the rock grains ( $var element = document.getElementById("moose-equation-6310a737-bd5e-4fd9-9da1-d5f1b7e4cc67");katex.render("c_{R}\\rho_{R}", element, {displayMode:false,throwOnError:false});$ ) are chosen to be constant.

Perfect, capped, weak-plane plasticity is used, so that the admissible zone is defined by $var element = document.getElementById("moose-equation-14f312aa-b768-4ff6-bb5f-b7c16b03bc64");katex.render("\\begin{array}{rcl} \\sigma_{zz} & \\leq & S_{T} \\ , \\\\ \\sigma_{zz} & \\geq & -S_{C} \\, \\\\ \\sqrt{\\sigma_{zx}^{2} + \\sigma_{zy}^{2}} + \\sigma_{zz}\\tan\\Phi & \\leq & C \\ . \\end{array}", element, {displayMode:true,throwOnError:false});$ Here $var element = document.getElementById("moose-equation-24102046-3109-4834-96e1-8ba011386074");katex.render("S_{T}", element, {displayMode:false,throwOnError:false});$ is the tensile strength, $var element = document.getElementById("moose-equation-75779aba-e282-4dee-aa23-40bc1fde3bf0");katex.render("S_{C}", element, {displayMode:false,throwOnError:false});$ is the compressive strength, $var element = document.getElementById("moose-equation-8292422a-f85f-4f98-9c8f-8133a45c5f1f");katex.render("C", element, {displayMode:false,throwOnError:false});$ is the cohesion and $var element = document.getElementById("moose-equation-393e4177-b4b6-4b65-aac5-aa434eb72011");katex.render("\\Phi", element, {displayMode:false,throwOnError:false});$ is the friction angle. The elastic tensor is chosen to be isotropic $var element = document.getElementById("moose-equation-f543a0f7-14f7-455c-ab0b-5f566e3c35b2");katex.render("E_{ijkl} = \\lambda \\delta_{ij}\\delta_{kl} + \\mu (\\delta_{ik}\\delta_{jl} + \\delta_{il}\\delta_{jk}) \\ .", element, {displayMode:true,throwOnError:false});$ The parameters in these expressions are chosen to be: $var element = document.getElementById("moose-equation-4310ae4c-2ece-4acb-ae79-7dca5de2f513");katex.render("S_{T}=1", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-e6e68c98-805e-4583-9061-dc8f0d46ec94");katex.render("S_{C}=1", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-9a42b5f2-b545-4152-9182-eef446c10f4b");katex.render("C=1", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-1e4e09f2-a954-4ae2-bed3-a6555d886c1e");katex.render("\\Phi=\\pi/4", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-019cc64a-6282-4bae-833e-7fd0d94557e8");katex.render("\\lambda=1/2", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-513cb2c1-1e5f-47b7-8105-29b164fd92e0");katex.render("\\mu=1/4", element, {displayMode:false,throwOnError:false});$ (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: $var element = document.getElementById("moose-equation-6790e229-a140-4ce1-89ee-247ce2a48759");katex.render("u_{z} = t \\ \\ \\ \\mathrm{for }\\ \\ z=1\\ ,", element, {displayMode:true,throwOnError:false});$ while the displacement on the element's bottom, and in the $var element = document.getElementById("moose-equation-a7535d63-5e38-47eb-9fd3-81f8e2d8fe92");katex.render("x", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-39add836-3866-423b-a3da-4bef89ad68b5");katex.render("y", element, {displayMode:false,throwOnError:false});$ 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<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 3
  xmin = -10
  xmax = 10
  zmin = 0
  zmax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y disp_z'
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [temperature]
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative<<<{"description": "Derivative of heat-energy-density wrt time", "href": "../../../../source/kernels/PorousFlowEnergyTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    base_name<<<{"description": "For mechanically-coupled systems, this Kernel will depend on the volumetric strain.  base_name should almost always be the same base_name as given to the TensorMechanics object that computes strain.  Supplying a base_name to this Kernel but not defining an associated TensorMechanics strain calculator means that this Kernel will not depend on volumetric strain.  That could be useful when models contain solid mechanics that is not coupled to porous flow, for example"}>>> = non_existent
  []
  [phe]
    type = PorousFlowPlasticHeatEnergy<<<{"description": "Plastic heat energy density source = (1 - porosity) * coeff * stress * plastic_strain_rate", "href": "../../../../source/kernels/PorousFlowPlasticHeatEnergy.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxKernels<<<{"href": "../../../../syntax/AuxKernels/index.html"}>>>]
  [disp_z]
    type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../../../source/auxkernels/FunctionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = disp_z
    function<<<{"description": "The function to use as the value"}>>> = z*t
  []
[]

[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."}>>> = temperature
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 0
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 0
  []
  [coh]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 100
  []
  [tanphi]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 1.0
  []
  [t_strength]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 1
  []
  [c_strength]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 1
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy<<<{"description": "This Material calculates the internal energy of solid rock grains, which is specific_heat_capacity * density * temperature.  Kernels multiply this by (1 - porosity) to find the energy density of the porous rock in a rock-fluid system", "href": "../../../../source/materials/PorousFlowMatrixInternalEnergy.html"}>>>
    specific_heat_capacity<<<{"description": "Specific heat capacity of the rock grains (J/kg/K)."}>>> = 2
    density<<<{"description": "Density of the rock grains"}>>> = 2
  []
  [temp]
    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"}>>>
    temperature<<<{"description": "Fluid temperature variable.  Note, the default is suitable if your simulation is using Kelvin units, but probably not for Celsius"}>>> = temperature
  []
  [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.2
  []
  [phe]
    type = ComputePlasticHeatEnergy<<<{"description": "Plastic heat energy density = stress * plastic_strain_rate", "href": "../../../../source/materials/ComputePlasticHeatEnergy.html"}>>>
  []
  [elasticity_tensor]
    type = ComputeElasticityTensor<<<{"description": "Compute an elasticity tensor.", "href": "../../../../source/materials/ComputeElasticityTensor.html"}>>>
    fill_method<<<{"description": "The fill method"}>>> = symmetric_isotropic
    C_ijkl<<<{"description": "Stiffness tensor for material"}>>> = '0.5 0.25'
  []
  [strain]
    type = ComputeIncrementalStrain<<<{"description": "Compute a strain increment and rotation increment for small strains.", "href": "../../../../source/materials/ComputeIncrementalStrain.html"}>>>
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
  [admissible]
    type = ComputeMultipleInelasticStress<<<{"description": "Compute state (stress and internal parameters such as plastic strains and internal parameters) using an iterative process.  Combinations of creep models and plastic models may be used.", "href": "../../../../source/materials/ComputeMultipleInelasticStress.html"}>>>
    inelastic_models<<<{"description": "The material objects to use to calculate stress and inelastic strains. Note: specify creep models first and plasticity models second."}>>> = mc
    perform_finite_strain_rotations<<<{"description": "Tensors are correctly rotated in finite-strain simulations.  For optimal performance you can set this to 'false' if you are only ever using small strains"}>>> = false
  []
  [mc]
    type = CappedWeakPlaneStressUpdate<<<{"description": "Capped weak-plane plasticity stress calculator", "href": "../../../../source/materials/CappedWeakPlaneStressUpdate.html"}>>>
    cohesion<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the cohesion.  Physically the cohesion should not be negative."}>>> = coh
    tan_friction_angle<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of tan(friction angle).  Physically the friction angle should be between 0 and 90deg."}>>> = tanphi
    tan_dilation_angle<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the tan(dilation angle).  Usually the dilation angle is not greater than the friction angle, and it is between 0 and 90deg."}>>> = tanphi
    tensile_strength<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the weak-plane tensile strength.  In physical situations this is positive (and always must be greater than negative compressive-strength."}>>> = t_strength
    compressive_strength<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the weak-plane compressive strength.  In physical situations this is positive."}>>> = c_strength
    tip_smoother<<<{"description": "The cone vertex at shear-stress = 0 will be smoothed by the given amount.  Typical value is 0.1*cohesion"}>>> = 0
    smoothing_tol<<<{"description": "Intersections of the yield surfaces will be smoothed by this amount (this is measured in units of stress).  Often this is related to other physical parameters (eg, 0.1*cohesion) but it is important to set this small enough so that the individual yield surfaces do not mix together in the smoothing process to produce a result where no stress is admissible (for example, mixing together tensile and compressive failure envelopes)."}>>> = 1
    yield_function_tol<<<{"description": "The return-map process will be deemed to have converged if all yield functions are within yield_function_tol of zero.  If this is set very low then precision-loss might be encountered: if the code detects precision loss then it also deems the return-map process has converged."}>>> = 1E-10
    perfect_guess<<<{"description": "Provide a guess to the Newton-Raphson procedure that is the result from perfect plasticity.  With severe hardening/softening this may be suboptimal."}>>> = true
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [temp]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = temperature
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [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\""}>>> = 'bcgs bjacobi 1E-15 1E-10 10000'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  dt = 1
  end_time = 10
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = tensile01
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

This implies the only non-zero component of total strain is $var element = document.getElementById("moose-equation-37318331-1fa0-4e34-ac08-68092a34101d");katex.render("\\epsilon^{\\mathrm{total}}_{zz} = t \\ .", element, {displayMode:true,throwOnError:false});$ The constitutive law implies $var element = document.getElementById("moose-equation-d556d94b-2bac-432d-bcb3-2153d7016c65");katex.render("\\sigma_{zz} = t \\ .", element, {displayMode:true,throwOnError:false});$ This stress is admissible for $var element = document.getElementById("moose-equation-b69c151d-2aaa-47a8-8af0-ab6469fcd9bb");katex.render("t\\leq 1", element, {displayMode:false,throwOnError:false});$ , while for $var element = document.getElementById("moose-equation-e3fd20fb-fc0c-4230-aba0-a53f44c872fe");katex.render("t>1", element, {displayMode:false,throwOnError:false});$ the system yields in tension: $var element = document.getElementById("moose-equation-744225f6-07d3-4222-affb-5c43459cb979");katex.render("\\sigma_{zz} = 1 \\ \\ \\ \\mathrm{for }\\ \\ t>1 \\ ,", element, {displayMode:true,throwOnError:false});$ and the plastic strain is, $var element = document.getElementById("moose-equation-ddf6214f-c315-49fb-a03b-28064bc5d019");katex.render("\\epsilon^{\\mathrm{plastic}}_{zz} = t - 1 \\ \\ \\ \\mathrm{for }\\ \\ t>1 \\ .", element, {displayMode:true,throwOnError:false});$ This means that the material's temperature should increase as $var element = document.getElementById("moose-equation-c17ccd88-5867-4502-853d-742870f79976");katex.render("c_{R}\\rho_{R}\\dot{T} = c \\ \\ \\ \\mathrm{for } t>1\\ \\ \\ ,", element, {displayMode:true,throwOnError:false});$ with no temperature increase for $var element = document.getElementById("moose-equation-8deaacf9-8814-43c9-b28a-adb414d15cae");katex.render("t\\leq 1", element, {displayMode:false,throwOnError:false});$ .

This result is reproduced exactly by PorousFlow.

Compressive failure

The top of the finite element is pushed downwards with displacement: $var element = document.getElementById("moose-equation-a2ef27ea-385c-40d8-90f1-6c2de99f99ff");katex.render("u_{z} = -t \\ \\ \\ \\mathrm{for }\\ \\ z=1\\ ,", element, {displayMode:true,throwOnError:false});$ while the displacement on the element's bottom, and in the $var element = document.getElementById("moose-equation-d299166a-8704-43d4-a6b1-3b21c63b00c2");katex.render("x", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-edf9bf47-45d4-49de-860f-d1b1b94881c9");katex.render("y", element, {displayMode:false,throwOnError:false});$ 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<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 3
  xmin = -10
  xmax = 10
  zmin = 0
  zmax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y disp_z'
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [temperature]
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative<<<{"description": "Derivative of heat-energy-density wrt time", "href": "../../../../source/kernels/PorousFlowEnergyTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    base_name<<<{"description": "For mechanically-coupled systems, this Kernel will depend on the volumetric strain.  base_name should almost always be the same base_name as given to the TensorMechanics object that computes strain.  Supplying a base_name to this Kernel but not defining an associated TensorMechanics strain calculator means that this Kernel will not depend on volumetric strain.  That could be useful when models contain solid mechanics that is not coupled to porous flow, for example"}>>> = non_existent
  []
  [phe]
    type = PorousFlowPlasticHeatEnergy<<<{"description": "Plastic heat energy density source = (1 - porosity) * coeff * stress * plastic_strain_rate", "href": "../../../../source/kernels/PorousFlowPlasticHeatEnergy.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxKernels<<<{"href": "../../../../syntax/AuxKernels/index.html"}>>>]
  [disp_z]
    type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../../../source/auxkernels/FunctionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = disp_z
    function<<<{"description": "The function to use as the value"}>>> = '-z*t'
  []
[]

[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."}>>> = temperature
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 0
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 0
  []
  [coh]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 100
  []
  [tanphi]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 1.0
  []
  [t_strength]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 1
  []
  [c_strength]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 1
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy<<<{"description": "This Material calculates the internal energy of solid rock grains, which is specific_heat_capacity * density * temperature.  Kernels multiply this by (1 - porosity) to find the energy density of the porous rock in a rock-fluid system", "href": "../../../../source/materials/PorousFlowMatrixInternalEnergy.html"}>>>
    specific_heat_capacity<<<{"description": "Specific heat capacity of the rock grains (J/kg/K)."}>>> = 2
    density<<<{"description": "Density of the rock grains"}>>> = 2
  []
  [temp]
    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"}>>>
    temperature<<<{"description": "Fluid temperature variable.  Note, the default is suitable if your simulation is using Kelvin units, but probably not for Celsius"}>>> = temperature
  []
  [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.2
  []
  [phe]
    type = ComputePlasticHeatEnergy<<<{"description": "Plastic heat energy density = stress * plastic_strain_rate", "href": "../../../../source/materials/ComputePlasticHeatEnergy.html"}>>>
  []
  [elasticity_tensor]
    type = ComputeElasticityTensor<<<{"description": "Compute an elasticity tensor.", "href": "../../../../source/materials/ComputeElasticityTensor.html"}>>>
    fill_method<<<{"description": "The fill method"}>>> = symmetric_isotropic
    C_ijkl<<<{"description": "Stiffness tensor for material"}>>> = '0.5 0.25'
  []
  [strain]
    type = ComputeIncrementalStrain<<<{"description": "Compute a strain increment and rotation increment for small strains.", "href": "../../../../source/materials/ComputeIncrementalStrain.html"}>>>
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
  [admissible]
    type = ComputeMultipleInelasticStress<<<{"description": "Compute state (stress and internal parameters such as plastic strains and internal parameters) using an iterative process.  Combinations of creep models and plastic models may be used.", "href": "../../../../source/materials/ComputeMultipleInelasticStress.html"}>>>
    inelastic_models<<<{"description": "The material objects to use to calculate stress and inelastic strains. Note: specify creep models first and plasticity models second."}>>> = mc
    perform_finite_strain_rotations<<<{"description": "Tensors are correctly rotated in finite-strain simulations.  For optimal performance you can set this to 'false' if you are only ever using small strains"}>>> = false
  []
  [mc]
    type = CappedWeakPlaneStressUpdate<<<{"description": "Capped weak-plane plasticity stress calculator", "href": "../../../../source/materials/CappedWeakPlaneStressUpdate.html"}>>>
    cohesion<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the cohesion.  Physically the cohesion should not be negative."}>>> = coh
    tan_friction_angle<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of tan(friction angle).  Physically the friction angle should be between 0 and 90deg."}>>> = tanphi
    tan_dilation_angle<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the tan(dilation angle).  Usually the dilation angle is not greater than the friction angle, and it is between 0 and 90deg."}>>> = tanphi
    tensile_strength<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the weak-plane tensile strength.  In physical situations this is positive (and always must be greater than negative compressive-strength."}>>> = t_strength
    compressive_strength<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the weak-plane compressive strength.  In physical situations this is positive."}>>> = c_strength
    tip_smoother<<<{"description": "The cone vertex at shear-stress = 0 will be smoothed by the given amount.  Typical value is 0.1*cohesion"}>>> = 0
    smoothing_tol<<<{"description": "Intersections of the yield surfaces will be smoothed by this amount (this is measured in units of stress).  Often this is related to other physical parameters (eg, 0.1*cohesion) but it is important to set this small enough so that the individual yield surfaces do not mix together in the smoothing process to produce a result where no stress is admissible (for example, mixing together tensile and compressive failure envelopes)."}>>> = 1
    yield_function_tol<<<{"description": "The return-map process will be deemed to have converged if all yield functions are within yield_function_tol of zero.  If this is set very low then precision-loss might be encountered: if the code detects precision loss then it also deems the return-map process has converged."}>>> = 1E-10
    perfect_guess<<<{"description": "Provide a guess to the Newton-Raphson procedure that is the result from perfect plasticity.  With severe hardening/softening this may be suboptimal."}>>> = true
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [temp]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = temperature
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [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\""}>>> = 'bcgs bjacobi 1E-15 1E-10 10000'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  dt = 1
  end_time = 10
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = compressive01
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

This implies only non-zero component of total strain is $var element = document.getElementById("moose-equation-67cfdf58-9147-4bf2-b369-d4e0752405a8");katex.render("\\epsilon^{\\mathrm{total}}_{zz} = -t \\ .", element, {displayMode:true,throwOnError:false});$ The constitutive law implies $var element = document.getElementById("moose-equation-62167f74-c87d-41ac-9a7b-589e9756ef5f");katex.render("\\sigma_{zz} = -t \\ .", element, {displayMode:true,throwOnError:false});$ This stress is admissible for $var element = document.getElementById("moose-equation-6ec47cae-3910-4d0c-86c0-4157e3779341");katex.render("t\\leq 1", element, {displayMode:false,throwOnError:false});$ , while for $var element = document.getElementById("moose-equation-4f814abb-590d-4b35-8613-d02884518522");katex.render("t>1", element, {displayMode:false,throwOnError:false});$ the system yields in compression $var element = document.getElementById("moose-equation-01cdaf58-ec9e-46da-8865-7b3545e05e03");katex.render("\\sigma_{zz} = -1 \\ \\ \\ \\mathrm{for }\\ \\ t>1 \\ ,", element, {displayMode:true,throwOnError:false});$ and the plastic strain is, $var element = document.getElementById("moose-equation-164d138c-384e-4ef2-83a1-d7ff3ee47913");katex.render("\\epsilon^{\\mathrm{plastic}}_{zz} = -(t - 1) \\ \\ \\ \\mathrm{for }\\ \\ t>1 \\ .", element, {displayMode:true,throwOnError:false});$ This means that the material's temperature should increase as $var element = document.getElementById("moose-equation-b71b32d9-dfdc-4b02-9dd4-523fb79e05b5");katex.render("c_{R}\\rho_{R}\\dot{T} = c \\ \\ \\ \\mathrm{for } t>1\\ \\ \\ ,", element, {displayMode:true,throwOnError:false});$ and there should be no temperature increase for $var element = document.getElementById("moose-equation-24eef9f4-fd10-4474-9899-3885e9897b40");katex.render("t\\leq 1", element, {displayMode:false,throwOnError:false});$ .

MOOSE produces this result exactly

Shear failure

The top of the finite element is sheared with displacement: $var element = document.getElementById("moose-equation-e105bdf7-2db4-48a5-9819-a7504d121759");katex.render("u_{x} = t \\ \\ \\ \\mathrm{for }\\ \\ z=1\\ ,", element, {displayMode:true,throwOnError:false});$ while the displacement on the element's bottom, and in the $var element = document.getElementById("moose-equation-87c1d217-a14c-4144-9eee-f73e25ff0edd");katex.render("y", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-98999309-c15d-4d93-bc3c-23186dee44bf");katex.render("z", element, {displayMode:false,throwOnError:false});$ 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<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 3
  xmin = -10
  xmax = 10
  zmin = 0
  zmax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y disp_z'
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [temperature]
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative<<<{"description": "Derivative of heat-energy-density wrt time", "href": "../../../../source/kernels/PorousFlowEnergyTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    base_name<<<{"description": "For mechanically-coupled systems, this Kernel will depend on the volumetric strain.  base_name should almost always be the same base_name as given to the TensorMechanics object that computes strain.  Supplying a base_name to this Kernel but not defining an associated TensorMechanics strain calculator means that this Kernel will not depend on volumetric strain.  That could be useful when models contain solid mechanics that is not coupled to porous flow, for example"}>>> = non_existent
  []
  [phe]
    type = PorousFlowPlasticHeatEnergy<<<{"description": "Plastic heat energy density source = (1 - porosity) * coeff * stress * plastic_strain_rate", "href": "../../../../source/kernels/PorousFlowPlasticHeatEnergy.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    coeff<<<{"description": "Heat energy density = coeff * stress * plastic_strain_rate"}>>> = 8
  []
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxKernels<<<{"href": "../../../../syntax/AuxKernels/index.html"}>>>]
  [disp_x]
    type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../../../source/auxkernels/FunctionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = disp_x
    function<<<{"description": "The function to use as the value"}>>> = 'z*t'
  []
[]

[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."}>>> = temperature
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 0
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 0
  []
  [coh]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 1
  []
  [tanphi]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 1.0
  []
  [tanpsi]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 0.0
  []
  [t_strength]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 10
  []
  [c_strength]
    type = TensorMechanicsHardeningConstant<<<{"description": "No hardening - the parameter is independent of the internal parameter(s)", "href": "../../../../source/userobjects/SolidMechanicsHardeningConstant.html"}>>>
    value<<<{"description": "The value of the parameter for all internal parameter.  This is perfect plasticity - there is no hardening."}>>> = 10
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy<<<{"description": "This Material calculates the internal energy of solid rock grains, which is specific_heat_capacity * density * temperature.  Kernels multiply this by (1 - porosity) to find the energy density of the porous rock in a rock-fluid system", "href": "../../../../source/materials/PorousFlowMatrixInternalEnergy.html"}>>>
    specific_heat_capacity<<<{"description": "Specific heat capacity of the rock grains (J/kg/K)."}>>> = 2
    density<<<{"description": "Density of the rock grains"}>>> = 2
  []
  [temp]
    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"}>>>
    temperature<<<{"description": "Fluid temperature variable.  Note, the default is suitable if your simulation is using Kelvin units, but probably not for Celsius"}>>> = temperature
  []
  [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.7
  []
  [phe]
    type = ComputePlasticHeatEnergy<<<{"description": "Plastic heat energy density = stress * plastic_strain_rate", "href": "../../../../source/materials/ComputePlasticHeatEnergy.html"}>>>
  []
  [elasticity_tensor]
    type = ComputeElasticityTensor<<<{"description": "Compute an elasticity tensor.", "href": "../../../../source/materials/ComputeElasticityTensor.html"}>>>
    fill_method<<<{"description": "The fill method"}>>> = symmetric_isotropic
    C_ijkl<<<{"description": "Stiffness tensor for material"}>>> = '0.5 0.25'
  []
  [strain]
    type = ComputeIncrementalStrain<<<{"description": "Compute a strain increment and rotation increment for small strains.", "href": "../../../../source/materials/ComputeIncrementalStrain.html"}>>>
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
  [admissible]
    type = ComputeMultipleInelasticStress<<<{"description": "Compute state (stress and internal parameters such as plastic strains and internal parameters) using an iterative process.  Combinations of creep models and plastic models may be used.", "href": "../../../../source/materials/ComputeMultipleInelasticStress.html"}>>>
    inelastic_models<<<{"description": "The material objects to use to calculate stress and inelastic strains. Note: specify creep models first and plasticity models second."}>>> = mc
    perform_finite_strain_rotations<<<{"description": "Tensors are correctly rotated in finite-strain simulations.  For optimal performance you can set this to 'false' if you are only ever using small strains"}>>> = false
  []
  [mc]
    type = CappedWeakPlaneStressUpdate<<<{"description": "Capped weak-plane plasticity stress calculator", "href": "../../../../source/materials/CappedWeakPlaneStressUpdate.html"}>>>
    cohesion<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the cohesion.  Physically the cohesion should not be negative."}>>> = coh
    tan_friction_angle<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of tan(friction angle).  Physically the friction angle should be between 0 and 90deg."}>>> = tanphi
    tan_dilation_angle<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the tan(dilation angle).  Usually the dilation angle is not greater than the friction angle, and it is between 0 and 90deg."}>>> = tanpsi
    tensile_strength<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the weak-plane tensile strength.  In physical situations this is positive (and always must be greater than negative compressive-strength."}>>> = t_strength
    compressive_strength<<<{"description": "A SolidMechanicsHardening UserObject that defines hardening of the weak-plane compressive strength.  In physical situations this is positive."}>>> = c_strength
    tip_smoother<<<{"description": "The cone vertex at shear-stress = 0 will be smoothed by the given amount.  Typical value is 0.1*cohesion"}>>> = 0
    smoothing_tol<<<{"description": "Intersections of the yield surfaces will be smoothed by this amount (this is measured in units of stress).  Often this is related to other physical parameters (eg, 0.1*cohesion) but it is important to set this small enough so that the individual yield surfaces do not mix together in the smoothing process to produce a result where no stress is admissible (for example, mixing together tensile and compressive failure envelopes)."}>>> = 1
    yield_function_tol<<<{"description": "The return-map process will be deemed to have converged if all yield functions are within yield_function_tol of zero.  If this is set very low then precision-loss might be encountered: if the code detects precision loss then it also deems the return-map process has converged."}>>> = 1E-10
    perfect_guess<<<{"description": "Provide a guess to the Newton-Raphson procedure that is the result from perfect plasticity.  With severe hardening/softening this may be suboptimal."}>>> = true
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [temp]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = temperature
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [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\""}>>> = 'bcgs bjacobi 1E-15 1E-10 10000'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  dt = 1
  end_time = 10
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = shear01
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

This implies only non-zero component of total strain is $var element = document.getElementById("moose-equation-f0200ae0-d114-4978-a6af-3b5375355616");katex.render("\\epsilon^{\\mathrm{total}}_{xz} = t \\ .", element, {displayMode:true,throwOnError:false});$ The constitutive law implies $var element = document.getElementById("moose-equation-2909bad5-c9b2-4df7-b0e0-2eace2245789");katex.render("\\sigma_{xz} = t/4 \\ .", element, {displayMode:true,throwOnError:false});$ This stress is admissible for $var element = document.getElementById("moose-equation-5b6f6b67-cb1c-4cc6-a526-23d66c540079");katex.render("t\\leq 4", element, {displayMode:false,throwOnError:false});$ , while for $var element = document.getElementById("moose-equation-e7754fb1-927b-4284-aa50-bcc3d09bfc14");katex.render("t>4", element, {displayMode:false,throwOnError:false});$ the system yields in shear: $var element = document.getElementById("moose-equation-0c83c9a7-9f97-4342-8d92-deb740014c10");katex.render("\\sigma_{xz} = 1 \\ \\ \\ \\mathrm{for }\\ \\ t>4 \\ ,", element, {displayMode:true,throwOnError:false});$ and the plastic strain is, $var element = document.getElementById("moose-equation-c66b42cf-e7b1-48c4-8c58-81b30631a19d");katex.render("\\epsilon^{\\mathrm{plastic}}_{xz} = t - 4 \\ \\ \\ \\mathrm{for }\\ \\ t>4 \\ .", element, {displayMode:true,throwOnError:false});$ This means that the material's temperature should increase as $var element = document.getElementById("moose-equation-5b031e3c-19c4-49cf-816a-e3ddce763427");katex.render("c_{R}\\rho_{R}\\dot{T} = c \\ \\ \\ \\mathrm{for }\\ \\ t>4 \\ ,", element, {displayMode:true,throwOnError:false});$ while the right-hand side is zero for $var element = document.getElementById("moose-equation-bed9d4e0-4819-4245-939d-c1a8265c0e56");katex.render("t\\leq 4", element, {displayMode:false,throwOnError:false});$ .

PorousFlow produces this result exactly.

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

# 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
    base_name = non_existent
  []
  [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 = ComputeIncrementalStrain
    displacements = 'disp_x disp_y disp_z'
  []
  [admissible]
    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
[]

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

# 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
    base_name = non_existent
  []
  [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 = ComputeIncrementalStrain
    displacements = 'disp_x disp_y disp_z'
  []
  [admissible]
    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
[]

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

# 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
    base_name = non_existent
  []
  [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 = ComputeIncrementalStrain
    displacements = 'disp_x disp_y disp_z'
  []
  [admissible]
    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
[]

Tensile failure
Compressive failure
Shear failure

Theory

Usage

Development

Heat energy produced by plastic deformation

Tensile failure

Compressive failure

Shear failure