Heat and fluid responses in 1D bars

Classic Newton cooling in a bar

Without fluids, mechanical deformation and sinks, the heat equation is $var element = document.getElementById("moose-equation-e7c49ac6-a6be-4cff-a0a9-38b5516663a1");katex.render("\\rho_{R}C_{R}(1-\\phi)\\dot{T} = \\nabla\\lambda\\nabla T \\ ,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-6c66c065-1e3c-4540-ada0-d61d8501942e");katex.render("\\phi", element, {displayMode:false,throwOnError:false});$ is the porosity, $var element = document.getElementById("moose-equation-fc675904-7133-40b9-bace-dcdd2f9164f1");katex.render("\\rho_{R}", element, {displayMode:false,throwOnError:false});$ is the rock grain density (kg.m $var element = document.getElementById("moose-equation-26662d13-6add-4ac3-81ea-b404720a050d");katex.render("^{-3}", element, {displayMode:false,throwOnError:false});$ ), $var element = document.getElementById("moose-equation-132bac96-d8a3-4f8d-a419-c20ff5a86a72");katex.render("C_{R}", element, {displayMode:false,throwOnError:false});$ is the rock grain specific heat capacity (J.kg $var element = document.getElementById("moose-equation-67bbd5b3-9f5c-489d-993a-9c128596c764");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ .K $var element = document.getElementById("moose-equation-bf0d465e-229c-4930-863b-f7785b1d4c54");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ), $var element = document.getElementById("moose-equation-3ad50910-0d03-4b7b-b527-ec800df04fde");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature, and $var element = document.getElementById("moose-equation-3ac2e8d9-9892-4b04-ab09-c271760b8ff3");katex.render("\\lambda", element, {displayMode:false,throwOnError:false});$ is the tensorial thermal conductivity of the porous material (J.s $var element = document.getElementById("moose-equation-940ff55f-caef-471c-b492-1cd1cf039978");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ .K $var element = document.getElementById("moose-equation-de972a01-0a47-4e1f-b47d-92e3c7b36367");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ .m $var element = document.getElementById("moose-equation-62313185-3f57-4825-997f-ac1286494b39");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ).

Below, the dynamics of this equation is explored, while this section concentrates on the steady-state situation. Consider the one-dimensional case where a bar sits between $var element = document.getElementById("moose-equation-08421d27-f387-4ab4-a18f-4758c1014873");katex.render("x=0", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-08b55681-a385-40c8-b950-541b80a5433c");katex.render("x=L", element, {displayMode:false,throwOnError:false});$ with a fixed temperature at $var element = document.getElementById("moose-equation-14af4f0a-8baa-4d80-8a94-5107756c02c8");katex.render("x=0", element, {displayMode:false,throwOnError:false});$ : $var element = document.getElementById("moose-equation-304db62e-b1e1-41f6-aa23-61156e86a121");katex.render("T(x=0, t) = T_{0} \\ ,", element, {displayMode:true,throwOnError:false});$ and a sink flux at the other end: $(1)var element = document.getElementById("moose-equation-a08cbf54-0d04-49e8-b983-36f0d764010d");katex.render("{\\mathrm{sink\\ strength}}\\, (J.m^{-2}.s^{-1}) = \\left.\\lambda\\frac{\\partial T}{\\partial x}\\right|_{x=L} = -C\\left(T - T_{e}\\right)_{x=L} \\ . ", element, {displayMode:true,throwOnError:false});$ Here $var element = document.getElementById("moose-equation-ffb03a79-5f80-40fe-ad63-904e55e0a884");katex.render("T_{e}", element, {displayMode:false,throwOnError:false});$ is a fixed quantity ("e" stands for "external"), and $var element = document.getElementById("moose-equation-97377bd4-fba0-482b-9caa-db856e6cd842");katex.render("C", element, {displayMode:false,throwOnError:false});$ is a constant conductance (J.m $var element = document.getElementById("moose-equation-45365387-a36f-47af-b991-d859768f043c");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ .s $var element = document.getElementById("moose-equation-b8549caa-cd8d-497f-a4db-34783ea1c2f6");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ .K $var element = document.getElementById("moose-equation-8ae844ae-53a3-4dc0-b290-6f7619b472a3");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ).

The solution is the linear function $var element = document.getElementById("moose-equation-ccebbc3f-e853-4305-8b13-a373c0daa14c");katex.render("T = T_{0} + \\frac{T_{e} - T_{0}}{\\lambda + C L} C x \\ .", element, {displayMode:true,throwOnError:false});$ The heat sink in Eq. (1) is a linear function of $var element = document.getElementById("moose-equation-716c91d9-2a78-4864-b5e8-467008ab0add");katex.render("T", element, {displayMode:false,throwOnError:false});$ , so the PorousFlowPiecewiseLinearSink may be employed.

The simulation is run in MOOSE using $var element = document.getElementById("moose-equation-ba293236-c9e4-4813-ab4d-eb0e8cf1e62d");katex.render("C=1", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-2e83d880-cad5-405b-a757-5a7654d388b8");katex.render("L=100", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-5d59f429-29bc-435a-99a4-e0250aa69867");katex.render("\\lambda=100", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-dc177149-e6b6-459e-8100-430d9eae487a");katex.render("T_{0}=2", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-c0f7ddb1-b317-4832-92a4-6f6990bbbc18");katex.render("T_{e}=1", element, {displayMode:false,throwOnError:false});$ . The solution is shown in Figure 1

Figure 1: The steady-state temperature in the bar. MOOSE agrees well with theory illustrating that piecewise-linear heat sinks/sources and heat conduction are correctly implemented in MOOSE.

# Newton cooling from a bar.  Heat conduction
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 100
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[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."}>>> = 'temp'
    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
  []
[]

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

[ICs<<<{"href": "../../../../syntax/ICs/index.html"}>>>]
  [temp]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = temp
    function<<<{"description": "The initial condition function."}>>> = '2-x/100'
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [conduction]
    type = PorousFlowHeatConduction<<<{"description": "Heat conduction in the Porous Flow module", "href": "../../../../source/kernels/PorousFlowHeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    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"}>>> = temp
  []
  [thermal_conductivity_irrelevant]
    type = PorousFlowThermalConductivityIdeal<<<{"description": "This Material calculates rock-fluid combined thermal conductivity by using a weighted sum.  Thermal conductivity = dry_thermal_conductivity + S^exponent * (wet_thermal_conductivity - dry_thermal_conductivity), where S is the aqueous saturation", "href": "../../../../source/materials/PorousFlowThermalConductivityIdeal.html"}>>>
    dry_thermal_conductivity<<<{"description": "The thermal conductivity of the rock matrix when the aqueous saturation is zero"}>>> = '1E2 0 0 0 1E2 0 0 0 1E2'
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [left]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 2
  []
  [newton]
    type = PorousFlowPiecewiseLinearSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.", "href": "../../../../source/bcs/PorousFlowPiecewiseLinearSink.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    pt_vals<<<{"description": "Tuple of pressure values (for the fluid_phase specified).  Must be monotonically increasing.  For heat fluxes that don't involve fluids, these are temperature values"}>>> = '0 1 2'
    multipliers<<<{"description": "Tuple of multiplying values.  The flux values are multiplied by these."}>>> = '-1 0 1'
    flux_function<<<{"description": "The flux.  The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE.  The functional form is useful for spatially or temporally varying sinks.  Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = 1
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [temp]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = temp
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0.5 0'
    end_point<<<{"description": "The ending of the line"}>>> = '100 0.5 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 11
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[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<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol '
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'gmres asm lu 100 NONZERO 2 1E-14 1E-12'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Steady
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = nc04
  execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = false
  [along_line]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    execute_vector_postprocessors_on<<<{"description": "Enable/disable the output of VectorPostprocessors"}>>> = timestep_end
  []
[]

Porepressure sink in a bar

These tests demonstrate that MOOSE behaves correctly when a simulation contains a sink. The sink is a piecewise linear function of pressure.

Darcy's equation for (single-phase) flow through a fully saturated medium without gravity and without sources is $var element = document.getElementById("moose-equation-71b86375-25a3-4968-824d-062599b3bca3");katex.render("\\frac{\\partial}{\\partial t}\\phi\\rho = \\nabla_{i}\\left(\\frac{\\rho \\kappa_{ij}}{\\mu} \\nabla_{j}P \\right) \\ ,", element, {displayMode:true,throwOnError:false});$ with the following notation:

$var element = document.getElementById("moose-equation-17c5f757-6fcd-4e4d-a5be-4af4a52e0a5f");katex.render("\\phi", element, {displayMode:false,throwOnError:false});$ is the medium's porosity;
$var element = document.getElementById("moose-equation-d15bf9f4-f655-4767-a712-30eb49091118");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ is the fluid density;
$var element = document.getElementById("moose-equation-b7bd8301-ad87-465a-938c-857a6af99684");katex.render("\\kappa_{ij}", element, {displayMode:false,throwOnError:false});$ is the permeability tensor;
$var element = document.getElementById("moose-equation-c0109d4b-0965-45b2-a7cc-3869e1fed906");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ is the fluid viscosity;
$var element = document.getElementById("moose-equation-ef6ec4c3-9248-40f5-a580-869fd26912d5");katex.render("\\partial/\\partial t", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-2143a3a0-e2a4-41b5-8c00-f28e37464e40");katex.render("\\nabla_{i}", element, {displayMode:false,throwOnError:false});$ denote the time and spatial derivatives, respectively.

Using $var element = document.getElementById("moose-equation-229f9962-f8cf-4075-b7cc-a6455c8917bf");katex.render("\\rho \\propto \\exp(P/B)", element, {displayMode:false,throwOnError:false});$ , where $var element = document.getElementById("moose-equation-267d2c3f-b442-473d-a997-59675d33e8ad");katex.render("B", element, {displayMode:false,throwOnError:false});$ is the fluid bulk modulus, Darcy's equation becomes $var element = document.getElementById("moose-equation-09670b58-5567-40fd-93ed-2a2edb24935f");katex.render("\\frac{\\partial}{\\partial t}\\rho = \\nabla_{i}\\alpha_{ij}\\nabla_{j}\\rho \\ ,", element, {displayMode:true,throwOnError:false});$ with $var element = document.getElementById("moose-equation-3d7833d5-d492-4e6f-bdab-09e0a987a0d6");katex.render("\\alpha_{ij} = \\frac{\\kappa_{ij}B}{\\mu\\phi} \\ .", element, {displayMode:true,throwOnError:false});$ Here the porosity and bulk modulus are assumed to be constant in space and time.

Consider the one-dimensional case where a bar sits between $var element = document.getElementById("moose-equation-d7b1941a-0834-4f74-863d-975b66f993dc");katex.render("x=0", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-04634428-9431-4bca-beab-2a77844bc78b");katex.render("x=L", element, {displayMode:false,throwOnError:false});$ with initial pressure distribution so $var element = document.getElementById("moose-equation-8bddb977-8d72-4e30-8faa-036844f3d962");katex.render("\\rho(x,t=0) = \\rho_{0}(x)", element, {displayMode:false,throwOnError:false});$ . Maintain the end $var element = document.getElementById("moose-equation-a12714d4-b72b-4c0f-bf95-1ad3e52bcee0");katex.render("x=0", element, {displayMode:false,throwOnError:false});$ at constant pressure, so that $var element = document.getElementById("moose-equation-e289967a-802b-45b2-a74a-18b29b5c493b");katex.render("\\rho(x=0, t) = \\rho_{0}(0)", element, {displayMode:false,throwOnError:false});$ . At the end $var element = document.getElementById("moose-equation-68afb9bc-4a43-4347-baf4-554b106d3f65");katex.render("x=L", element, {displayMode:false,throwOnError:false});$ , prescribe a sink flux $var element = document.getElementById("moose-equation-bc85b6de-5ec1-4311-bcbd-ad65cf4b2531");katex.render("\\left.\\frac{\\partial\\rho}{\\partial x}\\right|_{x=L} = -C\\left(\\rho - \\rho_{e}\\right)_{x=L} \\ ,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-490c7d42-b858-47f4-bd06-8461cb84927e");katex.render("\\rho_{e}", element, {displayMode:false,throwOnError:false});$ is a fixed quantity ("e" stands for "external"), and $var element = document.getElementById("moose-equation-a92d1ab4-6ffc-4cb4-b2d7-ca789eb5b4d3");katex.render("C", element, {displayMode:false,throwOnError:false});$ is a constant conductance. This corresponds to the flux $var element = document.getElementById("moose-equation-75cd13ef-d640-4f86-9e32-bc56e8bd607d");katex.render("\\left.\\frac{\\partial P}{\\partial x}\\right|_{x=L} = -CB\\left(1 - e^{(P_{e}-P)/B}\\right)_{x=L} \\ ,", element, {displayMode:true,throwOnError:false});$ which can easily be coded into a MOOSE input file: the flux is $var element = document.getElementById("moose-equation-3bf076ec-4303-46a9-b3a9-10839190fcdb");katex.render("\\rho\\kappa\\nabla P/\\mu = -CB\\kappa(e^{P/B} - e^{P_{e}/B})/\\mu", element, {displayMode:false,throwOnError:false});$ , and this may be represented by a piecewise linear function of pressure.

The solution of this problem is well known and is $var element = document.getElementById("moose-equation-d6b27812-31c2-48af-b7ca-d27d2d284875");katex.render("\\rho(x, t) = \\rho_{0}(0) - \\frac{\\rho_{0}(0) - \\rho_{e}}{1 + LC}Cx + \\sum_{n=1}^{\\infty} a_{n}\\sin \\frac{k_{n}x}{L}e^{-k_{n}^{2}\\alpha t/L^{2}} \\ ,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-5280598a-371d-4fe1-944a-698c2044f304");katex.render("k_{n}", element, {displayMode:false,throwOnError:false});$ is the $var element = document.getElementById("moose-equation-2a9334e0-b323-40c6-a98b-0f7ca5a26fa0");katex.render("n^{\\mathrm{th}}", element, {displayMode:false,throwOnError:false});$ positive root of the equation $var element = document.getElementById("moose-equation-cc4f830a-7694-4930-a3fa-9226d1c7ae9c");katex.render("LC\\tan k + k=0", element, {displayMode:false,throwOnError:false});$ ( $var element = document.getElementById("moose-equation-35b55529-3d90-4e77-9dd3-131e283d07f0");katex.render("k_{n}", element, {displayMode:false,throwOnError:false});$ is a little bigger than $var element = document.getElementById("moose-equation-8c6cc4c9-360b-45d5-b6bd-60e159f65b95");katex.render("(2n-1)\\pi/2", element, {displayMode:false,throwOnError:false});$ ), and $var element = document.getElementById("moose-equation-b199e42f-f268-4409-a92a-79d8a69425ab");katex.render("a_{n}", element, {displayMode:false,throwOnError:false});$ is determined from $var element = document.getElementById("moose-equation-46eae96b-8a5c-4198-94b7-405ab8c2ef3d");katex.render("a_{n}\\int_{0}^{L}\\sin^{2}\\frac{k_{n}x}{L}\\,\\mathrm{d}x = \\int_{0}^{L}\\left(\\rho_{0}(x) - \\rho_{0}(0) + \\frac{\\rho_{0}(0) - \\rho_{e}}{1 + LC}Cx\\right)\\sin \\frac{k_{n}x}{L}\\,\\mathrm{d}x \\ ,", element, {displayMode:true,throwOnError:false});$ which may be solved numerically (Mathematica is used to generate the solution in Figure 2).

The problem is solved in MOOSE using the following parameters:

Table 1: Parameter values used in the porepressure evolution tests

Parameter	Value
Bar length	$var element = document.getElementById("moose-equation-8c388b5c-ac98-4d41-b52d-a3eb0efbe3a6");katex.render("100\\,", element, {displayMode:false,throwOnError:false});$ m
Bar porosity	0.1
Bar permeability	$var element = document.getElementById("moose-equation-2e0a0eef-4629-49cc-ae15-0035db25efff");katex.render("10^{-15}\\,", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-d5a6843a-7208-436b-95df-655000a7667b");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$
Gravity	0
Water density	1000 $var element = document.getElementById("moose-equation-42c5b5fa-1b25-4e90-ba62-2b52e066b242");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ kg.m $var element = document.getElementById("moose-equation-ec86feca-420b-4664-8442-f7893776581f");katex.render("^{-3}", element, {displayMode:false,throwOnError:false});$
Water viscosity	0.001 $var element = document.getElementById("moose-equation-236e3b19-6719-440e-a81b-4623d2bab5b7");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ Pa.s
Water bulk modulus	1 $var element = document.getElementById("moose-equation-200b83ac-2863-4f06-b924-d7a764440a41");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ MPa
Initial porepressure $var element = document.getElementById("moose-equation-edafb031-b2c5-4e98-b484-b4ead0d73f55");katex.render("P_{0}", element, {displayMode:false,throwOnError:false});$	2 $var element = document.getElementById("moose-equation-0c52524f-89cb-4afb-912e-4ef3066e18f8");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ MPa
Environmental pressure $var element = document.getElementById("moose-equation-10254547-e2ea-4576-90ec-6336bed8b88a");katex.render("P_{e}", element, {displayMode:false,throwOnError:false});$	0
Conductance $var element = document.getElementById("moose-equation-19460ecb-c23a-4e12-b3b2-1442e2a3d355");katex.render("C", element, {displayMode:false,throwOnError:false});$	0.05389 $var element = document.getElementById("moose-equation-af000ec9-9aac-40cc-b407-73a4ffda1fb1");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-07b585e8-87bb-4dbe-af00-927b25b3c4db");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$

This conductance is chosen so at steadystate $var element = document.getElementById("moose-equation-7f9cb202-2924-4a53-ab03-c34373e07332");katex.render("\\rho(x=L)=2000\\,", element, {displayMode:false,throwOnError:false});$ kg.m $var element = document.getElementById("moose-equation-393aba01-61f0-4aa0-b2f8-312c6526fee8");katex.render("^{-3}", element, {displayMode:false,throwOnError:false});$ .

The problem is solved using 1000 elements along the $var element = document.getElementById("moose-equation-53611dd8-dc55-421a-a358-cfcbaa5cbc02");katex.render("x", element, {displayMode:false,throwOnError:false});$ direction ( $var element = document.getElementById("moose-equation-b7ca8713-8829-4a06-af2d-0e13a8261de3");katex.render("L=100\\,", element, {displayMode:false,throwOnError:false});$ m), and using 100 time-steps of size $var element = document.getElementById("moose-equation-0c98a450-9d99-46e2-b9b7-dc6d31175c76");katex.render("10^6\\,", element, {displayMode:false,throwOnError:false});$ s. Using fewer elements or fewer timesteps means the agreement with the theory is marginally poorer. Two tests are performed: one with transient flow, and one using the steadystate solver. In steady-state case the initial condition is $var element = document.getElementById("moose-equation-7e51cb06-0626-421c-93eb-166b5a98d16d");katex.render("P=2-x/L\\,", element, {displayMode:false,throwOnError:false});$ MPa, since the uniform $var element = document.getElementById("moose-equation-1746efe9-2ee0-4904-932b-518488525d8f");katex.render("P=2\\,", element, {displayMode:false,throwOnError:false});$ MPa does not converge. The results are shown in Figure 2.

Figure 2: The porepressure in the bar at 1E8s, and at steadystate. The pressure at is held fixed on the left-hand end, while the sink is applied at the right end. MOOSE agrees well with theory demonstrating that piecewise-linear sinks/sources and single-phase Darcy fluid flow are correctly implemented in MOOSE.

The transient simulation:

# Newton cooling from a bar.  1-phase transient
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 1000
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[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."}>>> = 'pressure'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.8
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-5
  []
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pressure]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 2E6
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
  []
  [flux]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 1e6
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    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"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pressure
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [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.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>> # irrelevant in this fully-saturated situation
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [left]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 2E6
  []
  [newton]
    type = PorousFlowPiecewiseLinearSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.", "href": "../../../../source/bcs/PorousFlowPiecewiseLinearSink.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    pt_vals<<<{"description": "Tuple of pressure values (for the fluid_phase specified).  Must be monotonically increasing.  For heat fluxes that don't involve fluids, these are temperature values"}>>> = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
    multipliers<<<{"description": "Tuple of multiplying values.  The flux values are multiplied by these."}>>> = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
    use_mobility<<<{"description": "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary.  In this case bare_flux is measured in Pa.m^-1.  This can be used in conjunction with other use_*"}>>> = false
    use_relperm<<<{"description": "If true, then fluxes are multiplied by relative permeability.  This can be used in conjunction with other use_*"}>>> = false
    fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy.  If not supplied, then this BC can only be a function of temperature"}>>> = 0
    flux_function<<<{"description": "The flux.  The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE.  The functional form is useful for spatially or temporally varying sinks.  Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = 1
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [porepressure]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = pressure
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0.5 0'
    end_point<<<{"description": "The ending of the line"}>>> = '100 0.5 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 20
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[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<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    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-12 1E-15 10000'
  []
[]

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

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = nc01
  [along_line]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    execute_vector_postprocessors_on<<<{"description": "Enable/disable the output of VectorPostprocessors"}>>> = final
  []
[]

The steady-state simulation:

# Newton cooling from a bar.  1-phase steady
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 1000
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[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."}>>> = 'pressure'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.8
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-5
  []
[]

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

[ICs<<<{"href": "../../../../syntax/ICs/index.html"}>>>]
  [pressure]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pressure
    function<<<{"description": "The initial condition function."}>>> = '(2-x/100)*1E6'
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [flux]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 1e6
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    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"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pressure
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [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.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>> # irrelevant in this fully-saturated situation
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [left]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 2E6
  []
  [newton]
    type = PorousFlowPiecewiseLinearSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.", "href": "../../../../source/bcs/PorousFlowPiecewiseLinearSink.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    pt_vals<<<{"description": "Tuple of pressure values (for the fluid_phase specified).  Must be monotonically increasing.  For heat fluxes that don't involve fluids, these are temperature values"}>>> = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
    multipliers<<<{"description": "Tuple of multiplying values.  The flux values are multiplied by these."}>>> = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
    use_mobility<<<{"description": "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary.  In this case bare_flux is measured in Pa.m^-1.  This can be used in conjunction with other use_*"}>>> = false
    use_relperm<<<{"description": "If true, then fluxes are multiplied by relative permeability.  This can be used in conjunction with other use_*"}>>> = false
    fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy.  If not supplied, then this BC can only be a function of temperature"}>>> = 0
    flux_function<<<{"description": "The flux.  The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE.  The functional form is useful for spatially or temporally varying sinks.  Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = 1
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [porepressure]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = pressure
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0.5 0'
    end_point<<<{"description": "The ending of the line"}>>> = '100 0.5 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 20
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  active<<<{"description": "If specified only the blocks named will be visited and made active"}>>> = 'andy'
  [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<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol '
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'gmres asm lu 100 NONZERO 2 1E-12 1E-15'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Steady
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = nc02
  execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = false
  [along_line]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    execute_vector_postprocessors_on<<<{"description": "Enable/disable the output of VectorPostprocessors"}>>> = timestep_end
  []
[]

Porepressure sink in a bar with heat

The simulation of the previous section is re-run, but this time heat flow is included. In this section it is assumed that the fluid specific enthalpy (J.kg $var element = document.getElementById("moose-equation-9c33a482-6ae4-4455-91a2-2da38d760d25");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ) is exactly equal to the fluid internal energy, and that internal energy is ideal: $var element = document.getElementById("moose-equation-9cc4d541-a42b-49f3-96ba-c4ccfc7f7ff4");katex.render("h = {\\mathcal E} = C_{v}T \\ .", element, {displayMode:true,throwOnError:false});$ This makes the arguments below simple without having to consider real fluids with complicated enthalpy and density expressions.

At the left end of the bar, the temperature is kept fixed: $var element = document.getElementById("moose-equation-60c0c447-89af-4a72-95c1-94a22cb4803a");katex.render("T(x=0, t) = T_{0} \\ .", element, {displayMode:true,throwOnError:false});$ At the other end of the bar, heat is removed only by the fluid flowing out of the system (see sink documentation). That is, there is a heat sink: $var element = document.getElementById("moose-equation-3064705a-541c-48b7-b53a-c81a132858e8");katex.render("{\\mathrm{sink\\ strength\\ }} (J.m^{-2}.s^{-1}) = -C {\\mathcal E}\\left(\\rho - \\rho_{e}\\right)_{x=L} \\ .", element, {displayMode:true,throwOnError:false});$ No other sinks or sources are applied to the heat equation.

With this setup, the steady-state temperature in the bar must be exactly $var element = document.getElementById("moose-equation-7fb0992e-917a-4055-9c2d-d790ccfde7fa");katex.render("T(x, t=\\infty) = T_{0} \\ .", element, {displayMode:true,throwOnError:false});$ For consider the fluid flowing from $var element = document.getElementById("moose-equation-5a73a76a-6bcc-4e54-84a6-b5e3bdd95be6");katex.render("x=0", element, {displayMode:false,throwOnError:false});$ to $var element = document.getElementById("moose-equation-aba4653c-db39-4437-9ecd-873078155cc9");katex.render("x=L", element, {displayMode:false,throwOnError:false});$ in order to assume steady-state. At $var element = document.getElementById("moose-equation-de0375b7-9de2-4823-8111-63403032a83c");katex.render("x=0", element, {displayMode:false,throwOnError:false});$ it must have temperature $var element = document.getElementById("moose-equation-265c9c0a-e1f4-47e2-a087-1be21e19fa34");katex.render("T_{0}", element, {displayMode:false,throwOnError:false});$ because that temperature is fixed at $var element = document.getElementById("moose-equation-5a39c149-4f62-478e-91ab-20681df754d3");katex.render("x=0", element, {displayMode:false,throwOnError:false});$ . It advects this temperature with it as it moves, so therefore at $var element = document.getElementById("moose-equation-275ed5d3-f919-427f-92b6-421c17cea42d");katex.render("t=\\infty", element, {displayMode:false,throwOnError:false});$ , this temperature has permeated throughout the entire bar. This occurs even without heat conduction, and is independent of the initial temperature of the bar.

MOOSE produces this result exactly. The input file is

# Newton cooling from a bar.  1-phase and heat, steady
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 100
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[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."}>>> = 'pressure temp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.8
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-5
  []
[]

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

[ICs<<<{"href": "../../../../syntax/ICs/index.html"}>>>]
  # have to start these reasonably close to their steady-state values
  [pressure]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pressure
    function<<<{"description": "The initial condition function."}>>> = '(2-x/100)*1E6'
  []
  [temperature]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = temp
    function<<<{"description": "The initial condition function."}>>> = 100+0.1*x
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [flux]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
  []
  [heat_advection]
    type = PorousFlowHeatAdvection<<<{"description": "Fully-upwinded heat flux, advected by the fluid", "href": "../../../../source/kernels/PorousFlowHeatAdvection.html"}>>>
    gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 1e6
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
    cv<<<{"description": "Constant specific heat capacity at constant volume (J/kg/K)"}>>> = 1e6
    porepressure_coefficient<<<{"description": "The enthalpy is internal_energy + P / density * porepressure_coefficient.  Physically this should be 1.0, but analytic solutions are simplified when it is zero"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    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"}>>> = temp
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pressure
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>> # irrelevant in this fully-saturated situation
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [leftp]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 2E6
  []
  [leftt]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 100
  []
  [newtonp]
    type = PorousFlowPiecewiseLinearSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.", "href": "../../../../source/bcs/PorousFlowPiecewiseLinearSink.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    pt_vals<<<{"description": "Tuple of pressure values (for the fluid_phase specified).  Must be monotonically increasing.  For heat fluxes that don't involve fluids, these are temperature values"}>>> = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
    multipliers<<<{"description": "Tuple of multiplying values.  The flux values are multiplied by these."}>>> = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
    use_mobility<<<{"description": "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary.  In this case bare_flux is measured in Pa.m^-1.  This can be used in conjunction with other use_*"}>>> = false
    use_relperm<<<{"description": "If true, then fluxes are multiplied by relative permeability.  This can be used in conjunction with other use_*"}>>> = false
    fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy.  If not supplied, then this BC can only be a function of temperature"}>>> = 0
    flux_function<<<{"description": "The flux.  The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE.  The functional form is useful for spatially or temporally varying sinks.  Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = 1
  []
  [newton]
    type = PorousFlowPiecewiseLinearSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.", "href": "../../../../source/bcs/PorousFlowPiecewiseLinearSink.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    pt_vals<<<{"description": "Tuple of pressure values (for the fluid_phase specified).  Must be monotonically increasing.  For heat fluxes that don't involve fluids, these are temperature values"}>>> = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
    multipliers<<<{"description": "Tuple of multiplying values.  The flux values are multiplied by these."}>>> = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
    use_mobility<<<{"description": "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary.  In this case bare_flux is measured in Pa.m^-1.  This can be used in conjunction with other use_*"}>>> = false
    use_relperm<<<{"description": "If true, then fluxes are multiplied by relative permeability.  This can be used in conjunction with other use_*"}>>> = false
    use_internal_energy<<<{"description": "If true, then fluxes are multiplied by fluid internal energy.  In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg).  This can be used in conjunction with other use_*"}>>> = true
    fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy.  If not supplied, then this BC can only be a function of temperature"}>>> = 0
    flux_function<<<{"description": "The flux.  The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE.  The functional form is useful for spatially or temporally varying sinks.  Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = 1
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [porepressure]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = pressure
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0.5 0'
    end_point<<<{"description": "The ending of the line"}>>> = '100 0.5 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 11
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [temperature]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = temp
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0.5 0'
    end_point<<<{"description": "The ending of the line"}>>> = '100 0.5 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 11
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[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<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol '
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'gmres asm lu 100 NONZERO 2 1E-8 1E-15'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Steady
  solve_type = Newton
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = nc06
  execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  [along_line]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    execute_vector_postprocessors_on<<<{"description": "Enable/disable the output of VectorPostprocessors"}>>> = timestep_end
  []
[]

Hot ideal fluid in a bar

This test uses a similar setup to the previous section, except that here an ideal fluid is used. The use of an ideal gas simplifies the equations. Only the steady-state is studied in this section.

The governing equation for the fluid's porepressure $var element = document.getElementById("moose-equation-2e7c90a4-bf52-4168-8f11-b865643daa5c");katex.render("P", element, {displayMode:false,throwOnError:false});$ is $var element = document.getElementById("moose-equation-043f94da-1209-49d8-bec5-58c9c4376180");katex.render("\\nabla \\frac{\\rho\\kappa}{\\mu}\\nabla P = 0 \\ .", element, {displayMode:true,throwOnError:false});$ It is assumed that $var element = document.getElementById("moose-equation-77682ffa-1727-4d3f-ba9f-2dbe95c4bb07");katex.render("\\kappa", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-2684f7fa-e4af-40c1-9051-e07841eca5fc");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ are constant, and that $var element = document.getElementById("moose-equation-68956073-57c5-4ef5-91b4-c9558a3358d6");katex.render("\\rho = \\frac{MP}{RT} \\ ,", element, {displayMode:true,throwOnError:false});$ holds (this is the ideal gas equation of state). In this formula $var element = document.getElementById("moose-equation-72603649-0e8b-45d9-bc88-ea256684dc9d");katex.render("M", element, {displayMode:false,throwOnError:false});$ is the gas molar mass, $var element = document.getElementById("moose-equation-da91ed17-7b95-4cf7-bf56-4201b0c655d4");katex.render("R", element, {displayMode:false,throwOnError:false});$ is the gas constant and $var element = document.getElementById("moose-equation-26d6fedc-1976-4955-9faf-a0836aa4f3c7");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature.

The equation governing the temperature is assumed to be just the fluid advection equation $var element = document.getElementById("moose-equation-bd1ddb18-91da-438e-97bb-7ba4ca29f9d3");katex.render("\\nabla\\frac{h\\rho\\kappa}{\\mu}\\nabla P = 0 \\ .", element, {displayMode:true,throwOnError:false});$ As in the previous section, heat conduction could be added, but it is actually irrelevant since the solution to the problem below is constant $var element = document.getElementById("moose-equation-318b505c-d544-4eb5-8368-6d4d18e4add6");katex.render("T", element, {displayMode:false,throwOnError:false});$ . The enthalpy, $var element = document.getElementById("moose-equation-587bcaa2-9a11-415f-b6a3-3a58442d634d");katex.render("h", element, {displayMode:false,throwOnError:false});$ , for an ideal gas is $var element = document.getElementById("moose-equation-78ed6f77-f047-4a92-a7e6-8e565d210da7");katex.render("h = C_{v}T + \\frac{P}{\\rho} = C_{v}T + \\frac{RT}{M} = C_{p}T \\ .", element, {displayMode:true,throwOnError:false});$

The boundary conditions at the left-hand end are $var element = document.getElementById("moose-equation-d0ce7dd3-6ceb-42d9-b3f7-eed0839206c7");katex.render("P(x=0) = P_{0} \\ \\ \\ \\mathrm{and}\\ \\ \\ T(x=0)=T_{0} \\ .", element, {displayMode:true,throwOnError:false});$ Physically these correspond to fluid and heat being removed or added to the left-hand end by some external source in order to keep the porepressure and temperature fixed.

The porepressure boundary condition at the right-hand end of the bar is $(2)var element = document.getElementById("moose-equation-fab53cff-3ccc-400f-bfdc-4a49fcddd387");katex.render("\\mathrm{sink\\ flux\\ } (kg.m^{-2}.s^{-1}) = \\left. \\frac{\\rho\\kappa}{\\mu}\\nabla P \\right|_{x=L} = \\left. -C\\frac{\\rho\\kappa}{\\mu} (P - P_{e}) \\right|_{x=L} \\ . ", element, {displayMode:true,throwOnError:false});$ Physically this corresponds to the mass-flow through the boundary being proportional to $var element = document.getElementById("moose-equation-2d5525bd-a884-4dce-a8a9-57ba42f40478");katex.render("P-P_{e}", element, {displayMode:false,throwOnError:false});$ . Here $var element = document.getElementById("moose-equation-109e139d-f143-41cd-b738-df8351786f08");katex.render("P_{e}", element, {displayMode:false,throwOnError:false});$ is a fixed "environmental" porepressure, and this acts as a source or sink of fluid. $var element = document.getElementById("moose-equation-25148109-2875-4b38-afa5-2d0659767e12");katex.render("C", element, {displayMode:false,throwOnError:false});$ is the "conductance" of the boundary. Notice the appearance of $var element = document.getElementById("moose-equation-a0da98f0-4d62-4fcd-acea-4f0f3e4307de");katex.render("\\rho \\kappa/\\mu", element, {displayMode:false,throwOnError:false});$ in the LHS of this equation means that this is truly a flux of fluid mass (measured in kg.m $var element = document.getElementById("moose-equation-ea293b85-ce05-4a7c-916b-91834ebf26fb");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ .s $var element = document.getElementById("moose-equation-aaa43e2b-c0c9-486d-8724-c923da890d0d");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ), and the appearance of $var element = document.getElementById("moose-equation-e3b8d123-adb8-42c2-b14e-e62d8f0f11eb");katex.render("\\rho\\kappa/\\mu", element, {displayMode:false,throwOnError:false});$ on the RHS means that a PorousFlowPiecewiseLinearFlux may be used with use_mobility=true.

The temperature boundary condition at the right-hand end of the bar is $var element = document.getElementById("moose-equation-773abc40-5123-43e4-b320-9de130df2c93");katex.render("\\mathrm{heat\\ flux\\ } (J.m^{-2}.s^{-1}) = \\left. \\frac{h\\rho\\kappa}{\\mu}\\nabla P \\right|_{x=L} = \\left. -C\\frac{h\\rho\\kappa}{\\mu} (P - P_{e}) \\right|_{x=L} \\ .", element, {displayMode:true,throwOnError:false});$ Comparing this with Eq. (2), it is seen that this is exactly the heat loss (or gain) at the boundary corresponding to the loss (or gain) of the fluid. Notice the appearance of $var element = document.getElementById("moose-equation-d94140eb-fd4f-4868-9894-b29ac39b3973");katex.render("h\\rho \\kappa/\\mu", element, {displayMode:false,throwOnError:false});$ in the LHS of this equation means that this is truly a flux of fluid mass (measured in J.m $var element = document.getElementById("moose-equation-d522e97c-5c26-4677-96ba-6be55a7b4f0a");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ .s $var element = document.getElementById("moose-equation-42f3469e-26a7-4703-9533-62e4f5ec89f6");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ), and the appearance of $var element = document.getElementById("moose-equation-ebfbee68-b723-4131-95cb-101209d323ae");katex.render("h\\rho\\kappa/\\mu", element, {displayMode:false,throwOnError:false});$ on the RHS means that a PorousFlowPiecewiseLinearFlux may be used with use_mobility=true and use_enthalpy=true.

There is a clear similarity between the fluid and heat equations. The heat equation does not actually depend on temperature, and is simply $var element = document.getElementById("moose-equation-7313779d-05fa-443f-b806-3db27984d86b");katex.render("0 = \\nabla (P\\nabla P) \\ ,", element, {displayMode:true,throwOnError:false});$ which is solved by $var element = document.getElementById("moose-equation-22521fe1-c5f3-49aa-84d8-6bdce254f140");katex.render("P^{2} = P_{0}^{2} + Ax \\ .", element, {displayMode:true,throwOnError:false});$ The fluid equation then yields $var element = document.getElementById("moose-equation-0a169b6a-4373-43ca-b1ee-e0c29723a5db");katex.render("T(x) = T_{0} \\ .", element, {displayMode:true,throwOnError:false});$ The constant $var element = document.getElementById("moose-equation-8eb809c3-fc05-4503-8eaa-bc8a5fc6fc1b");katex.render("A", element, {displayMode:false,throwOnError:false});$ may be determined from the either of the boundary conditions. For the special case of $var element = document.getElementById("moose-equation-94e8b981-0cf7-471b-aeed-238929c2d989");katex.render("P_{e}=0", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-a9e0fd18-f277-4fb0-9c37-3ed131b3a55a");katex.render("2LC=1", element, {displayMode:false,throwOnError:false});$ , the solution is $var element = document.getElementById("moose-equation-c47f708e-a003-4ac8-8ec4-78823ff71e28");katex.render("P = P_{0} \\sqrt{1 - \\frac{x}{2L}} \\ .", element, {displayMode:true,throwOnError:false});$ MOOSE produces this result exactly, as illustrated in Figure 3

Figure 3: The steady-state porepressure and temperature distributions in the bar. MOOSE agrees well with theory illustrating that piecewise-linear fluid and heat sinks/sources as well as ideal fluids are correctly implemented in MOOSE.

# Newton cooling from a bar.  1-phase ideal fluid and heat, steady
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 100
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[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."}>>> = 'pressure temp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.8
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-5
  []
[]

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

[ICs<<<{"href": "../../../../syntax/ICs/index.html"}>>>]
  # have to start these reasonably close to their steady-state values
  [pressure]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pressure
    function<<<{"description": "The initial condition function."}>>> = '200-0.5*x'
  []
  [temperature]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = temp
    function<<<{"description": "The initial condition function."}>>> = 180+0.1*x
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [flux]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
  []
  [heat_advection]
    type = PorousFlowHeatAdvection<<<{"description": "Fully-upwinded heat flux, advected by the fluid", "href": "../../../../source/kernels/PorousFlowHeatAdvection.html"}>>>
    gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [idealgas]
    type = IdealGasFluidProperties<<<{"description": "Fluid properties for an ideal gas", "href": "../../../../source/fluidproperties/IdealGasFluidProperties.html"}>>>
    molar_mass<<<{"description": "Constant molar mass of the fluid (kg/mol)"}>>> = 1.4
    gamma<<<{"description": "gamma value (cp/cv)"}>>> = 1.2
    mu<<<{"description": "Dynamic viscosity, Pa.s"}>>> = 1.2
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    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"}>>> = temp
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pressure
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [dens0]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = idealgas
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1.1 0 0 0 1.1 0 0 0 1.1'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>> # irrelevant in this fully-saturated situation
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [leftp]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 200
  []
  [leftt]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 180
  []
  [newtonp]
    type = PorousFlowPiecewiseLinearSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.", "href": "../../../../source/bcs/PorousFlowPiecewiseLinearSink.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    pt_vals<<<{"description": "Tuple of pressure values (for the fluid_phase specified).  Must be monotonically increasing.  For heat fluxes that don't involve fluids, these are temperature values"}>>> = '-200 0 200'
    multipliers<<<{"description": "Tuple of multiplying values.  The flux values are multiplied by these."}>>> = '-200 0 200'
    use_mobility<<<{"description": "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary.  In this case bare_flux is measured in Pa.m^-1.  This can be used in conjunction with other use_*"}>>> = true
    use_relperm<<<{"description": "If true, then fluxes are multiplied by relative permeability.  This can be used in conjunction with other use_*"}>>> = true
    fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy.  If not supplied, then this BC can only be a function of temperature"}>>> = 0
    flux_function<<<{"description": "The flux.  The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE.  The functional form is useful for spatially or temporally varying sinks.  Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = 0.005 # 1/2/L
  []
  [newtont]
    type = PorousFlowPiecewiseLinearSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.", "href": "../../../../source/bcs/PorousFlowPiecewiseLinearSink.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    pt_vals<<<{"description": "Tuple of pressure values (for the fluid_phase specified).  Must be monotonically increasing.  For heat fluxes that don't involve fluids, these are temperature values"}>>> = '-200 0 200'
    multipliers<<<{"description": "Tuple of multiplying values.  The flux values are multiplied by these."}>>> = '-200 0 200'
    use_mobility<<<{"description": "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary.  In this case bare_flux is measured in Pa.m^-1.  This can be used in conjunction with other use_*"}>>> = true
    use_relperm<<<{"description": "If true, then fluxes are multiplied by relative permeability.  This can be used in conjunction with other use_*"}>>> = true
    use_enthalpy<<<{"description": "If true, then fluxes are multiplied by enthalpy.  In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg).  This can be used in conjunction with other use_*"}>>> = true
    fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy.  If not supplied, then this BC can only be a function of temperature"}>>> = 0
    flux_function<<<{"description": "The flux.  The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE.  The functional form is useful for spatially or temporally varying sinks.  Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = 0.005 # 1/2/L
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [porepressure]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = pressure
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0.5 0'
    end_point<<<{"description": "The ending of the line"}>>> = '100 0.5 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 11
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [temperature]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = temp
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0.5 0'
    end_point<<<{"description": "The ending of the line"}>>> = '100 0.5 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 11
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[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
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Steady
  solve_type = Newton
  nl_rel_tol = 1E-10
  nl_abs_tol = 1E-15
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = nc08
  execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  [along_line]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    execute_vector_postprocessors_on<<<{"description": "Enable/disable the output of VectorPostprocessors"}>>> = timestep_end
  []
[]

(moose/modules/porous_flow/test/tests/newton_cooling/nc04.i)

# Newton cooling from a bar.  Heat conduction
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 100
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp'
    number_fluid_phases = 0
    number_fluid_components = 0
  []
[]

[Variables]
  [temp]
  []
[]

[ICs]
  [temp]
    type = FunctionIC
    variable = temp
    function = '2-x/100'
  []
[]

[Kernels]
  [conduction]
    type = PorousFlowHeatConduction
    variable = temp
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [thermal_conductivity_irrelevant]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '1E2 0 0 0 1E2 0 0 0 1E2'
  []
[]

[BCs]
  [left]
    type = DirichletBC
    variable = temp
    boundary = left
    value = 2
  []
  [newton]
    type = PorousFlowPiecewiseLinearSink
    variable = temp
    boundary = right
    pt_vals = '0 1 2'
    multipliers = '-1 0 1'
    flux_function = 1
  []
[]

[VectorPostprocessors]
  [temp]
    type = LineValueSampler
    variable = temp
    start_point = '0 0.5 0'
    end_point = '100 0.5 0'
    sort_by = x
    num_points = 11
    execute_on = timestep_end
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol '
    petsc_options_value = 'gmres asm lu 100 NONZERO 2 1E-14 1E-12'
  []
[]

[Executioner]
  type = Steady
[]

[Outputs]
  file_base = nc04
  execute_on = timestep_end
  exodus = false
  [along_line]
    type = CSV
    execute_vector_postprocessors_on = timestep_end
  []
[]

(moose/modules/porous_flow/test/tests/newton_cooling/nc01.i)

# Newton cooling from a bar.  1-phase transient
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 1000
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pressure'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.8
    alpha = 1e-5
  []
[]

[Variables]
  [pressure]
    initial_condition = 2E6
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pressure
  []
  [flux]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    gravity = '0 0 0'
    variable = pressure
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 1e6
    density0 = 1000
    thermal_expansion = 0
    viscosity = 1e-3
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey # irrelevant in this fully-saturated situation
    n = 2
    phase = 0
  []
[]

[BCs]
  [left]
    type = DirichletBC
    variable = pressure
    boundary = left
    value = 2E6
  []
  [newton]
    type = PorousFlowPiecewiseLinearSink
    variable = pressure
    boundary = right
    pt_vals = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
    multipliers = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 1
  []
[]

[VectorPostprocessors]
  [porepressure]
    type = LineValueSampler
    variable = pressure
    start_point = '0 0.5 0'
    end_point = '100 0.5 0'
    sort_by = x
    num_points = 20
    execute_on = timestep_end
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E8
  dt = 1E6
[]

[Outputs]
  file_base = nc01
  [along_line]
    type = CSV
    execute_vector_postprocessors_on = final
  []
[]

(moose/modules/porous_flow/test/tests/newton_cooling/nc02.i)

# Newton cooling from a bar.  1-phase steady
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 1000
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pressure'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.8
    alpha = 1e-5
  []
[]

[Variables]
  [pressure]
  []
[]

[ICs]
  [pressure]
    type = FunctionIC
    variable = pressure
    function = '(2-x/100)*1E6'
  []
[]

[Kernels]
  [flux]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    gravity = '0 0 0'
    variable = pressure
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 1e6
    density0 = 1000
    thermal_expansion = 0
    viscosity = 1e-3
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey # irrelevant in this fully-saturated situation
    n = 2
    phase = 0
  []
[]

[BCs]
  [left]
    type = DirichletBC
    variable = pressure
    boundary = left
    value = 2E6
  []
  [newton]
    type = PorousFlowPiecewiseLinearSink
    variable = pressure
    boundary = right
    pt_vals = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
    multipliers = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 1
  []
[]

[VectorPostprocessors]
  [porepressure]
    type = LineValueSampler
    variable = pressure
    start_point = '0 0.5 0'
    end_point = '100 0.5 0'
    sort_by = x
    num_points = 20
    execute_on = timestep_end
  []
[]

[Preconditioning]
  active = 'andy'
  [andy]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol '
    petsc_options_value = 'gmres asm lu 100 NONZERO 2 1E-12 1E-15'
  []
[]

[Executioner]
  type = Steady
[]

[Outputs]
  file_base = nc02
  execute_on = timestep_end
  exodus = false
  [along_line]
    type = CSV
    execute_vector_postprocessors_on = timestep_end
  []
[]

(moose/modules/porous_flow/test/tests/newton_cooling/nc06.i)

# Newton cooling from a bar.  1-phase and heat, steady
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 100
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pressure temp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.8
    alpha = 1e-5
  []
[]

[Variables]
  [pressure]
  []
  [temp]
  []
[]

[ICs]
  # have to start these reasonably close to their steady-state values
  [pressure]
    type = FunctionIC
    variable = pressure
    function = '(2-x/100)*1E6'
  []
  [temperature]
    type = FunctionIC
    variable = temp
    function = 100+0.1*x
  []
[]

[Kernels]
  [flux]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    gravity = '0 0 0'
    variable = pressure
  []
  [heat_advection]
    type = PorousFlowHeatAdvection
    gravity = '0 0 0'
    variable = temp
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 1e6
    density0 = 1000
    thermal_expansion = 0
    viscosity = 1e-3
    cv = 1e6
    porepressure_coefficient = 0
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey # irrelevant in this fully-saturated situation
    n = 2
    phase = 0
  []
[]

[BCs]
  [leftp]
    type = DirichletBC
    variable = pressure
    boundary = left
    value = 2E6
  []
  [leftt]
    type = DirichletBC
    variable = temp
    boundary = left
    value = 100
  []
  [newtonp]
    type = PorousFlowPiecewiseLinearSink
    variable = pressure
    boundary = right
    pt_vals = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
    multipliers = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 1
  []
  [newton]
    type = PorousFlowPiecewiseLinearSink
    variable = temp
    boundary = right
    pt_vals = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
    multipliers = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
    use_mobility = false
    use_relperm = false
    use_internal_energy = true
    fluid_phase = 0
    flux_function = 1
  []
[]

[VectorPostprocessors]
  [porepressure]
    type = LineValueSampler
    variable = pressure
    start_point = '0 0.5 0'
    end_point = '100 0.5 0'
    sort_by = x
    num_points = 11
    execute_on = timestep_end
  []
  [temperature]
    type = LineValueSampler
    variable = temp
    start_point = '0 0.5 0'
    end_point = '100 0.5 0'
    sort_by = x
    num_points = 11
    execute_on = timestep_end
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol '
    petsc_options_value = 'gmres asm lu 100 NONZERO 2 1E-8 1E-15'
  []
[]

[Executioner]
  type = Steady
  solve_type = Newton
[]

[Outputs]
  file_base = nc06
  execute_on = timestep_end
  [along_line]
    type = CSV
    execute_vector_postprocessors_on = timestep_end
  []
[]

(moose/modules/porous_flow/test/tests/newton_cooling/nc08.i)

# Newton cooling from a bar.  1-phase ideal fluid and heat, steady
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 100
  ny = 1
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pressure temp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.8
    alpha = 1e-5
  []
[]

[Variables]
  [pressure]
  []
  [temp]
  []
[]

[ICs]
  # have to start these reasonably close to their steady-state values
  [pressure]
    type = FunctionIC
    variable = pressure
    function = '200-0.5*x'
  []
  [temperature]
    type = FunctionIC
    variable = temp
    function = 180+0.1*x
  []
[]

[Kernels]
  [flux]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    gravity = '0 0 0'
    variable = pressure
  []
  [heat_advection]
    type = PorousFlowHeatAdvection
    gravity = '0 0 0'
    variable = temp
  []
[]

[FluidProperties]
  [idealgas]
    type = IdealGasFluidProperties
    molar_mass = 1.4
    gamma = 1.2
    mu = 1.2
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [dens0]
    type = PorousFlowSingleComponentFluid
    fp = idealgas
    phase = 0
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.1 0 0 0 1.1 0 0 0 1.1'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey # irrelevant in this fully-saturated situation
    n = 2
    phase = 0
  []
[]

[BCs]
  [leftp]
    type = DirichletBC
    variable = pressure
    boundary = left
    value = 200
  []
  [leftt]
    type = DirichletBC
    variable = temp
    boundary = left
    value = 180
  []
  [newtonp]
    type = PorousFlowPiecewiseLinearSink
    variable = pressure
    boundary = right
    pt_vals = '-200 0 200'
    multipliers = '-200 0 200'
    use_mobility = true
    use_relperm = true
    fluid_phase = 0
    flux_function = 0.005 # 1/2/L
  []
  [newtont]
    type = PorousFlowPiecewiseLinearSink
    variable = temp
    boundary = right
    pt_vals = '-200 0 200'
    multipliers = '-200 0 200'
    use_mobility = true
    use_relperm = true
    use_enthalpy = true
    fluid_phase = 0
    flux_function = 0.005 # 1/2/L
  []
[]

[VectorPostprocessors]
  [porepressure]
    type = LineValueSampler
    variable = pressure
    start_point = '0 0.5 0'
    end_point = '100 0.5 0'
    sort_by = x
    num_points = 11
    execute_on = timestep_end
  []
  [temperature]
    type = LineValueSampler
    variable = temp
    start_point = '0 0.5 0'
    end_point = '100 0.5 0'
    sort_by = x
    num_points = 11
    execute_on = timestep_end
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Steady
  solve_type = Newton
  nl_rel_tol = 1E-10
  nl_abs_tol = 1E-15
[]

[Outputs]
  file_base = nc08
  execute_on = timestep_end
  [along_line]
    type = CSV
    execute_vector_postprocessors_on = timestep_end
  []
[]

Classic Newton cooling in a bar
Porepressure sink in a bar
Porepressure sink in a bar with heat
Hot ideal fluid in a bar

Theory

Usage

Development

Heat and fluid responses in 1D bars

Classic Newton cooling in a bar

Porepressure sink in a bar

Porepressure sink in a bar with heat

Hot ideal fluid in a bar