Step 1 - First Thermal Mortar Contact
Continuing from the final step in the mechanical contact tutorial, for the first time we incorporate thermal contact, staying within the mortar method.
#
# A first attempt at thermo mechanical contact
# https://mooseframework.inl.gov/modules/combined/tutorials/introduction/step01.html
#
[GlobalParams]
displacements = 'disp_x disp_y'
block = 0
[]
[Mesh]
[generated1]
type = GeneratedMeshGenerator
dim = 2
nx = 5
ny = 15
xmin = -0.6
xmax = -0.1
ymax = 5
bias_y = 0.9
boundary_name_prefix = pillar1
[]
[generated2]
type = GeneratedMeshGenerator
dim = 2
nx = 6
ny = 15
xmin = 0.1
xmax = 0.6
ymax = 4.999
bias_y = 0.9
boundary_name_prefix = pillar2
boundary_id_offset = 4
[]
[collect_meshes]
type = MeshCollectionGenerator
inputs = 'generated1 generated2'
[]
patch_update_strategy = iteration
[]
[Variables]
# temperature field variable
[T]
# initialize to an average temperature
initial_condition = 50
order = FIRST
family = LAGRANGE
[]
# temperature lagrange multiplier
[Tlm]
block = 'pillars_secondary_subdomain'
order = FIRST
family = LAGRANGE
[]
[]
[Kernels]
[heat_conduction]
type = HeatConduction
variable = T
[]
[dTdt]
type = HeatConductionTimeDerivative
variable = T
[]
[]
[Modules/TensorMechanics/Master]
[all]
add_variables = true
strain = FINITE
generate_output = 'vonmises_stress'
[]
[]
[Contact]
[pillars]
primary = pillar1_right
secondary = pillar2_left
model = frictionless
formulation = mortar
[]
[]
[Constraints]
# thermal contact constraint
[Tlm]
type = GapConductanceConstraint
variable = Tlm
secondary_variable = T
use_displaced_mesh = true
k = 1e-1
primary_boundary = pillar1_right
primary_subdomain = pillars_primary_subdomain
secondary_boundary = pillar2_left
secondary_subdomain = pillars_secondary_subdomain
[]
[]
[BCs]
[bottom_x]
type = DirichletBC
variable = disp_x
boundary = 'pillar1_bottom pillar2_bottom'
value = 0
[]
[bottom_y]
type = DirichletBC
variable = disp_y
boundary = 'pillar1_bottom pillar2_bottom'
value = 0
[]
[Pressure]
[sides]
boundary = 'pillar1_left pillar2_right'
function = 1e4*t^2
[]
[]
# thermal boundary conditions (pillars are heated/cooled from the bottom)
[heat_left]
type = DirichletBC
variable = T
boundary = pillar1_bottom
value = 100
[]
[cool_right]
type = DirichletBC
variable = T
boundary = pillar2_bottom
value = 0
[]
[]
[Materials]
[elasticity]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1e9
poissons_ratio = 0.3
[]
[stress]
type = ComputeFiniteStrainElasticStress
[]
# thermal properties
[thermal_conductivity]
type = HeatConductionMaterial
thermal_conductivity = 100
specific_heat = 1
[]
[density]
type = Density
density = 1
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
line_search = none
# we deal with the saddle point structure of the system by adding a small shift
petsc_options_iname = '-pc_type -pc_factor_shift_type'
petsc_options_value = 'lu nonzero'
end_time = 5
dt = 0.1
[Predictor]
type = SimplePredictor
scale = 1
[]
[]
[Outputs]
exodus = true
print_linear_residuals = false
perf_graph = true
[]
Note that we do not have a "thermal contact action", so we need to add all objects required for mortar based thermal contact manually.
Input file
Variables
For the thermal part of the problem we add the temperature T variable and supply an initial_condition of 50 temperature units (this value is right between the 0 and 100 values we're setting below for the cooling and heating BCs). Note that like the other physics variables and objects this variable is restricted to block 0 (see GlobalParams).
The second variable we add is the Lagrange multiplier for the thermal contact. Its physical meaning is the heat flux! This variable will be block restricted to the lower dimensional mortar subdomain the [Contact] action created for us (the name of that subdomain is the name of the contact action subblock with _secondary_subdomain appended).
The "family" and "order" parameters define the shape function family and element order. First order Lagrange (no connection to "Lagrange multiplier!") are the MOOSE default settings. In future inputs you might see those parameters omitted. Lagrange shape functions are purely nodal interpolary shape functions.
Constraints
The GapConductanceConstraint is the bare bones gap conductance constraint object available in MOOSE. It could serve as a starting point for the development of more complex gap conductance models (taking into account surface roughness and contact pressure for example).
The constraint object acts on the Lagrange multiplier "variable" Tlm, and it affects the "secondary_variable" T, the temperature, enforcing the energy balance.
We make sure to enable the evaluation of the constraint on the displaced mesh to capture large displacement and deformatio by setting the "use_displaced_mesh" parameter to true. The displaced mesh is a copy of the mesh that has all nodes moved by their displacement variable values. Try setting this parameter to false and compare the results (hint the gap distance the constraint sees will be computed from the original mesh and will not change over time).
We chose a gap conductance "k". The heat flux across the gap will be computed as where is the width of the gap.
Then we set the primary and secondary boundaries and subdomains. Again for the subdomains we choose the lower dimensional meshes that the mechanical contact action set up for us. If you have a purely thermal contact problem you will need to create those lower dimensional subdomains manually using the LowerDBlockFromSidesetGenerator mesh generator.
BCs
heat_left and cool_right are two Dirichlet boundary conditions that set the bottom of the left cantilever to an elevated temperature and the bottom of the right cantilever to a lowered temperature, so we can establish a temperature difference and a resulting heat flux between the two cantilevers.
Materials
We need to add two material blocks for thermal properties. HeatConductionMaterial provides thermal conductivity and specific heat. And the Density material computes an updated material density on the displaced mesh, taking volumetric deformation into account. The specific_heat and density properties are used by the HeatConductionTimeDerivative kernel, and the thermal_conductivity property is used by the HeatConduction kernel.
Executioner
We make a small change to the PETSc options and add the -pc_factor_shift_type
nonzero option in "petsc_options_iname" and "petsc_options_value" respectively. This helps the direct solver deal with the saddle point structure of the Jacobian matrix (with zeroes on the diagonal).
Note that we also lower the timestep to better see what is going on. You can set it back to 0.5 to convince yourself that the problem still converges just as well.
Once you've answered the questions and run this example we will move on to Step 2 in which we introduce a third block participating in the contact problem.