Coupled Thermomechanics Bar Problem

Overview

This problem is one of 1D, coupled heat transfer and solid mechanics. The domain is a bar fixed at the left end with the right end free to displace due to thermal expansion. The left end of the bar has a set temperature, and the right end has a flux boundary condition where the flux depends on the total displacement of the bar. Thus, the displacement depends on temperature, which in turn depends on displacement.

Bar Problem Details

$var element = document.getElementById("moose-equation-3dbdbab9-3ffd-42c3-8f63-a7e4ec5129c3");katex.render("x", element, {displayMode:false,throwOnError:false});$ is the x-coordinate
$var element = document.getElementById("moose-equation-5d6b359a-f366-4d3d-9965-84868c7650ca");katex.render("u", element, {displayMode:false,throwOnError:false});$ is displacement
$var element = document.getElementById("moose-equation-80668f90-a4cf-4ee2-bdb8-242f35d66b77");katex.render("L", element, {displayMode:false,throwOnError:false});$ is the length (10)
$var element = document.getElementById("moose-equation-1bc70600-e59c-47fe-b6e7-268592477e5c");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$ is the coefficient of thermal expansion (1E-6)
$var element = document.getElementById("moose-equation-45a5fba0-73ee-44f4-8d21-18b6a4ec49a5");katex.render("T", element, {displayMode:false,throwOnError:false});$ is temperature
$var element = document.getElementById("moose-equation-02d679e1-7d62-490d-8c34-ca9bd9025d10");katex.render("T_0", element, {displayMode:false,throwOnError:false});$ is the initial temperature of the bar (300)
$var element = document.getElementById("moose-equation-b81ab32c-045a-4484-b4e8-fd3b667debf3");katex.render("T(0)", element, {displayMode:false,throwOnError:false});$ is the fixed temperature at the left end of the bar ( $var element = document.getElementById("moose-equation-30d13824-6a5f-4688-be62-cee9748ac99e");katex.render("T_0", element, {displayMode:false,throwOnError:false});$ )
$var element = document.getElementById("moose-equation-40efc44d-b15d-44c2-b0d5-aec4b2d3ecd7");katex.render("k", element, {displayMode:false,throwOnError:false});$ is the thermal conductivity (1000)
$var element = document.getElementById("moose-equation-d02ebf4e-ec8c-4fd2-a70f-5b7b62649dbd");katex.render("\\beta", element, {displayMode:false,throwOnError:false});$ is a constant (28500)
$var element = document.getElementById("moose-equation-49746c8e-2130-4228-866d-11c5de504b2a");katex.render("\\gamma", element, {displayMode:false,throwOnError:false});$ is a constant (1E6)
$var element = document.getElementById("moose-equation-2b18e1e0-91b1-4b9a-9285-7542e2b06eab");katex.render("q", element, {displayMode:false,throwOnError:false});$ is the flux ( $var element = document.getElementById("moose-equation-76345057-4dfb-4b28-b685-d2d5b4653498");katex.render("-\\beta-\\gamma u(L)", element, {displayMode:false,throwOnError:false});$ )

With a fixed temperature at one end and a flux boundary condition at the other, the temperature in the bar will vary linearly. The displacement at the end of the bar is then $var element = document.getElementById("moose-equation-dcafdcf5-e61f-4184-a2a4-f53104498ab2");katex.render("\\begin{aligned} u(L) =& \\int_0^L\\alpha(T(x)-T_0) \\mathrm{d}x\\\\ =& \\int_0^L\\alpha(T(0)+\\frac{T(L)-T(0)}{L}x-T_0) \\mathrm{d}x\\\\ =& \\int_0^L\\alpha(T_0-T_0+\\frac{T(L)-T(0)}{L}x) \\mathrm{d}x\\\\ =& \\int_0^L\\alpha(\\frac{T(L)-T(0)}{L}x) \\mathrm{d}x\\\\ =& \\frac{\\alpha}{2L}(T(L)-T(0))x^2\\rvert_0^L\\\\ =& \\frac{\\alpha L}{2}(T(L)-T(0)) \\end{aligned}", element, {displayMode:true,throwOnError:false});$

The differential equation for heat transfer in the bar is $var element = document.getElementById("moose-equation-df63841c-2e91-47b6-bbe0-6d3dccbc684d");katex.render("kT''=0", element, {displayMode:true,throwOnError:false});$ which after integration becomes $var element = document.getElementById("moose-equation-2b1c7534-6fe9-4967-9cdc-60acc3d9cf6b");katex.render("kT + C_1x + C_2 = 0", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-4bb84206-e8c1-443a-bc49-51479430802b");katex.render("C_1", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-69a4f9e8-da1a-426f-84e8-3ce10134d142");katex.render("C_2", element, {displayMode:false,throwOnError:false});$ are constants of integration.

$var element = document.getElementById("moose-equation-dbefc502-6cf5-4a84-8647-17d16c61f945");katex.render("T(0)", element, {displayMode:false,throwOnError:false});$ = $var element = document.getElementById("moose-equation-5cdd886b-9438-4015-a38b-869fb495e049");katex.render("T_0", element, {displayMode:false,throwOnError:false});$ , which leads to $var element = document.getElementById("moose-equation-aa248897-7d73-4c55-9d85-deffc2a8c4f9");katex.render("C_2=-kT_0", element, {displayMode:false,throwOnError:false});$ .

We also have $var element = document.getElementById("moose-equation-ca7c538c-ebce-45e9-926d-84340022fd42");katex.render("-kT=C_1=q", element, {displayMode:false,throwOnError:false});$ .

To solve for T(L), we have $var element = document.getElementById("moose-equation-1b1d4170-ad7a-4425-9da5-5eca7500b16b");katex.render("\\begin{aligned} T(L) =& T_0+\\left(\\frac{\\beta L}{k} + \\frac{\\gamma\\alpha L^2}{2k}(T(L)-T(0))\\right)\\\\ T(L) =& \\frac{\\gamma\\alpha L^2}{2k}T(L) + \\frac{\\beta L}{k} + T_0 - \\frac{\\gamma\\alpha L^2}{2k}T_0\\\\ T(L) =& \\frac{\\frac{\\beta L}{k}+\\left(1-\\frac{\\gamma\\alpha L^2}{2k}\\right)T_0}{1-\\frac{\\gamma\\alpha L^2}{2k}} \\end{aligned}", element, {displayMode:true,throwOnError:false});$

Table 1: Results for bar with coupled thermomechanics.

	Analytic Value	Bison Value	Error
u(L)	1.5E-3	1.5E-3	0%
q	-3E4	-3E4	0%
T(L)	600	600	0%

Input Syntax

#
# Quantity  Analytical
# --------  ----------
#   u(L)      1.5e-3
#   q           -3e4
#   T(L)         600
#
[Mesh<<<{"href": "../../../syntax/Mesh/index.html"}>>>]
  use_displaced_mesh<<<{"description": "Create the displaced mesh if the 'displacements' parameter is set. If this is 'false', a displaced mesh will not be created, regardless of whether 'displacements' is set."}>>> = false
  [bar]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 3
    nx<<<{"description": "Number of elements in the X direction"}>>> = 10
    ny<<<{"description": "Number of elements in the Y direction"}>>> = 1
    nz<<<{"description": "Number of elements in the Z direction"}>>> = 1
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 10
  []
[]

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

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [unity]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1
  []
[]

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

[Kernels<<<{"href": "../../../syntax/Kernels/index.html"}>>>]
  [heat]
    type = HeatConduction<<<{"description": "Diffusive heat conduction term $-\\nabla\\cdot(k\\nabla T)$ of the thermal energy conservation equation", "href": "../../../source/kernels/HeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
  []
[]

[Physics<<<{"href": "../../../syntax/Physics/index.html"}>>>]
  [SolidMechanics<<<{"href": "../../../syntax/Physics/SolidMechanics/index.html"}>>>]
    [QuasiStatic<<<{"href": "../../../syntax/Physics/SolidMechanics/QuasiStatic/index.html"}>>>]
      [all]
        strain<<<{"description": "Strain formulation"}>>> = SMALL
        incremental<<<{"description": "Use incremental or total strain (if not explicitly specified this defaults to incremental for finite strain and total for small strain)"}>>> = true
        add_variables<<<{"description": "Add the displacement variables"}>>> = true
        eigenstrain_names<<<{"description": "List of eigenstrains to be applied in this strain calculation"}>>> = eigenstrain
        generate_output<<<{"description": "Add scalar quantity output for stress and/or strain"}>>> = 'strain_xx strain_yy strain_zz'
        temperature<<<{"description": "The temperature"}>>> = temperature
      []
    []
  []
[]

[BCs<<<{"href": "../../../syntax/NuclearMaterials/BCs/index.html"}>>>]
  [x]
    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"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [y]
    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"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'top bottom'
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [z_bot]
    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"}>>> = disp_z
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'top bottom'
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [T]
    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"}>>> = temperature
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left'
    value<<<{"description": "Value of the BC"}>>> = 300
  []
  [flux]
    type = MatNeumannBC<<<{"description": "Imposes the integrated boundary condition $\\frac{C \\partial u}{\\partial n}=M*h$, where $h$ is a constant, $M$ is a material property, and $C$ is a coefficient defined by the kernel for $u$.", "href": "../../../source/bcs/MatNeumannBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    value<<<{"description": "For a Laplacian problem, the value of the gradient dotted with the normals on the boundary."}>>> = 1
    boundary_material<<<{"description": "Material property multiplying the constant that will be enforced by the BC"}>>> = flux_material
  []
[]

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../../../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 1
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0
  []
  [small_stress]
    type = ComputeFiniteStrainElasticStress<<<{"description": "Compute stress using elasticity for finite strains", "href": "../../../source/materials/ComputeFiniteStrainElasticStress.html"}>>>
  []
  [thermal_expansion_strain]
    type = ComputeThermalExpansionEigenstrain<<<{"description": "Computes eigenstrain due to thermal expansion with a constant coefficient", "href": "../../../source/materials/ComputeThermalExpansionEigenstrain.html"}>>>
    stress_free_temperature<<<{"description": "Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation"}>>> = 300
    thermal_expansion_coeff<<<{"description": "Thermal expansion coefficient"}>>> = 1e-6
    temperature<<<{"description": "Coupled temperature"}>>> = temperature
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = eigenstrain
  []
  [k]
    type = ParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = thermal_conductivity
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 1000
  []
  [flux]
    type = ParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = flux_material
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = disp_x
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '28500 + 1e6 * disp_x'
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = right
  []
[]

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

  l_max_its = 50
  nl_max_its = 50
  nl_rel_tol = 1e-14
  nl_abs_tol = 1e-10
  l_tol = 1e-9

  start_time = 0.0
  end_time = 1.0
  dt = 1.0
  dtmin = 1.0
[]

[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [analytical_T]
    type = ParsedPostprocessor<<<{"description": "Computes a parsed expression with post-processors", "href": "../../../source/postprocessors/ParsedPostprocessor.html"}>>>
    pp_names<<<{"description": "Post-processors arguments"}>>> = ''
    constant_names<<<{"description": "Vector of constants used in the parsed function (use this for kB etc.)"}>>> = 'beta L k gamma alpha T0'
    constant_expressions<<<{"description": "Vector of values for the constants in constant_names (can be an FParser expression)"}>>> = '28500 10 1000 1e6 1e-6 300'
    expression<<<{"description": "function expression"}>>> = 'gammaTerm := 1-gamma*alpha*L*L/2/k;
                (beta*L/k + gammaTerm*T0)/gammaTerm'
  []
  [analytical_u]
    type = ParsedPostprocessor<<<{"description": "Computes a parsed expression with post-processors", "href": "../../../source/postprocessors/ParsedPostprocessor.html"}>>>
    pp_names<<<{"description": "Post-processors arguments"}>>> = 'analytical_T'
    constant_names<<<{"description": "Vector of constants used in the parsed function (use this for kB etc.)"}>>> = 'L alpha T0'
    constant_expressions<<<{"description": "Vector of values for the constants in constant_names (can be an FParser expression)"}>>> = '10 1e-6 300'
    expression<<<{"description": "function expression"}>>> = 'alpha*L/2*(analytical_T-T0)'
  []
  [analytical_flux]
    type = ParsedPostprocessor<<<{"description": "Computes a parsed expression with post-processors", "href": "../../../source/postprocessors/ParsedPostprocessor.html"}>>>
    pp_names<<<{"description": "Post-processors arguments"}>>> = 'analytical_u'
    constant_names<<<{"description": "Vector of constants used in the parsed function (use this for kB etc.)"}>>> = 'beta gamma'
    constant_expressions<<<{"description": "Vector of values for the constants in constant_names (can be an FParser expression)"}>>> = '28500 1e6'
    expression<<<{"description": "function expression"}>>> = '-beta-gamma*analytical_u'
  []
  [bison_u]
    type = SideAverageValue<<<{"description": "Computes the average value of a variable on a sideset. Note that this cannot be used on the centerline of an axisymmetric model.", "href": "../../../source/postprocessors/SideAverageValue.html"}>>>
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
  []
  [bison_T]
    type = SideAverageValue<<<{"description": "Computes the average value of a variable on a sideset. Note that this cannot be used on the centerline of an axisymmetric model.", "href": "../../../source/postprocessors/SideAverageValue.html"}>>>
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = temperature
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
  []
  [bison_flux]
    type = SideDiffusiveFluxAverage<<<{"description": "Computes the integral of the diffusive flux over the specified boundary", "href": "../../../source/postprocessors/SideDiffusiveFluxAverage.html"}>>>
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = temperature
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    diffusivity<<<{"description": "The name of the diffusivity material property that will be used in the flux computation. This must be provided if the variable is of finite element type"}>>> = thermal_conductivity
  []
  [error_u]
    type = ParsedPostprocessor<<<{"description": "Computes a parsed expression with post-processors", "href": "../../../source/postprocessors/ParsedPostprocessor.html"}>>>
    expression<<<{"description": "function expression"}>>> = '(bison_u-analytical_u)/analytical_u*100'
    pp_names<<<{"description": "Post-processors arguments"}>>> = 'bison_u analytical_u'
  []
  [error_T]
    type = ParsedPostprocessor<<<{"description": "Computes a parsed expression with post-processors", "href": "../../../source/postprocessors/ParsedPostprocessor.html"}>>>
    expression<<<{"description": "function expression"}>>> = '(bison_T-analytical_T)/analytical_T*100'
    pp_names<<<{"description": "Post-processors arguments"}>>> = 'bison_T analytical_T'
  []
  [error_flux]
    type = ParsedPostprocessor<<<{"description": "Computes a parsed expression with post-processors", "href": "../../../source/postprocessors/ParsedPostprocessor.html"}>>>
    expression<<<{"description": "function expression"}>>> = '(bison_flux-analytical_flux)/analytical_flux*100'
    pp_names<<<{"description": "Post-processors arguments"}>>> = 'bison_flux analytical_flux'
  []
[]

(assessment/verification/Thermomechanics/bar/bar.i)

#
# Quantity  Analytical
# --------  ----------
#   u(L)      1.5e-3
#   q           -3e4
#   T(L)         600
#
[Mesh]
  use_displaced_mesh = false
  [bar]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 10
    ny = 1
    nz = 1
    xmax = 10
  []
[]

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

[AuxVariables]
  [unity]
    initial_condition = 1
  []
[]

[Variables]
  [temperature]
    initial_condition = 300
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
    use_displaced_mesh = false
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        incremental = true
        add_variables = true
        eigenstrain_names = eigenstrain
        generate_output = 'strain_xx strain_yy strain_zz'
        temperature = temperature
      []
    []
  []
[]

[BCs]
  [x]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  []
  [y]
    type = DirichletBC
    variable = disp_y
    boundary = 'top bottom'
    value = 0.0
  []
  [z_bot]
    type = DirichletBC
    variable = disp_z
    boundary = 'top bottom'
    value = 0.0
  []
  [T]
    type = DirichletBC
    variable = temperature
    boundary = 'left'
    value = 300
  []
  [flux]
    type = MatNeumannBC
    variable = temperature
    boundary = right
    value = 1
    boundary_material = flux_material
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0
  []
  [small_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [thermal_expansion_strain]
    type = ComputeThermalExpansionEigenstrain
    stress_free_temperature = 300
    thermal_expansion_coeff = 1e-6
    temperature = temperature
    eigenstrain_name = eigenstrain
  []
  [k]
    type = ParsedMaterial
    property_name = thermal_conductivity
    expression = 1000
  []
  [flux]
    type = ParsedMaterial
    property_name = flux_material
    coupled_variables = disp_x
    expression = '28500 + 1e6 * disp_x'
    boundary = right
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton

  l_max_its = 50
  nl_max_its = 50
  nl_rel_tol = 1e-14
  nl_abs_tol = 1e-10
  l_tol = 1e-9

  start_time = 0.0
  end_time = 1.0
  dt = 1.0
  dtmin = 1.0
[]

[Outputs]
  csv = true
[]

[Postprocessors]
  [analytical_T]
    type = ParsedPostprocessor
    pp_names = ''
    constant_names = 'beta L k gamma alpha T0'
    constant_expressions = '28500 10 1000 1e6 1e-6 300'
    expression = 'gammaTerm := 1-gamma*alpha*L*L/2/k;
                (beta*L/k + gammaTerm*T0)/gammaTerm'
  []
  [analytical_u]
    type = ParsedPostprocessor
    pp_names = 'analytical_T'
    constant_names = 'L alpha T0'
    constant_expressions = '10 1e-6 300'
    expression = 'alpha*L/2*(analytical_T-T0)'
  []
  [analytical_flux]
    type = ParsedPostprocessor
    pp_names = 'analytical_u'
    constant_names = 'beta gamma'
    constant_expressions = '28500 1e6'
    expression = '-beta-gamma*analytical_u'
  []
  [bison_u]
    type = SideAverageValue
    variable = disp_x
    boundary = right
  []
  [bison_T]
    type = SideAverageValue
    variable = temperature
    boundary = right
  []
  [bison_flux]
    type = SideDiffusiveFluxAverage
    variable = temperature
    boundary = right
    diffusivity = thermal_conductivity
  []
  [error_u]
    type = ParsedPostprocessor
    expression = '(bison_u-analytical_u)/analytical_u*100'
    pp_names = 'bison_u analytical_u'
  []
  [error_T]
    type = ParsedPostprocessor
    expression = '(bison_T-analytical_T)/analytical_T*100'
    pp_names = 'bison_T analytical_T'
  []
  [error_flux]
    type = ParsedPostprocessor
    expression = '(bison_flux-analytical_flux)/analytical_flux*100'
    pp_names = 'bison_flux analytical_flux'
  []
[]

Overview
Input Syntax