Solution without any constraint along Z-direction

[chakra34]

Hi,

I am having a weird observation and was hoping to get some insight from the community. I am showing a simplified version of the problem as an example here.

So, I have a quarter cylinder geometry for which I am solving the thermo-elastic problem under a certain temperature field. The boundary conditions are such that the face normal to global X-axis has Ux=0, Y-axis has Uy=0, and Z-axis has Uz=0. This is well constrained problem and the solution looks like below:

The corresponding moose_input file looks like so:

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = input_geometry.e
  []
[]
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Functions]
  [./temp_function]
    type  = ParsedFunction
    value = '200.0*(z+0)'
  [../]
[]
[BCs]
  [face_fix_y]
    type = DirichletBC
    variable = disp_y
    boundary = fix_y
    value = 0
  []
  [face_fix_x]
    type = DirichletBC
    variable = disp_x
    boundary = fix_x
    value = 0
  []
  [face_fix_z]
    type = DirichletBC
    variable = disp_z
    boundary = fix_z
    value = 0
  []
[]
[Modules]
  [TensorMechanics]
    [Master/All]
      strain = small
      incremental = False
      add_variables = True
      temperature = temp
      eigenstrain_names = eigenstrain
      generate_output = 'stress_xx stress_yy stress_zz'
      use_automatic_differentiation = False
    []
  []
[]
[Materials]
  [elasticity_tensor_core]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 150.0e9
    poissons_ratio = 0.3
    block = 1
  []
  [elasticity_tensor_skin]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 150.0e9
    poissons_ratio = 0.3
    block = 2
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [thermal_expansion_strain_core]
    type = ComputeThermalExpansionEigenstrain
    thermal_expansion_coeff = 1.0e-4
    stress_free_temperature = 100.0
    temperature = temp
    eigenstrain_name = eigenstrain
    block = 1
  []
  [thermal_expansion_strain_skin]
    type = ComputeThermalExpansionEigenstrain
    thermal_expansion_coeff = 1.0e-4
    stress_free_temperature = 100.0
    temperature = temp
    eigenstrain_name = eigenstrain
    block = 2
  []
[]
[AuxVariables]
  [temp]
    order = FIRST
    family = LAGRANGE
  []
  [stress_rr_core]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_tt_core]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zz_core]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_rr_skin]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_tt_skin]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zz_skin]
    order = CONSTANT
    family = MONOMIAL
  []
  [eigenstrain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [eigenstrain_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [eigenstrain_zz]
    order = CONSTANT
    family = MONOMIAL
  []
[]
[AuxKernels]
  [set_temp]
    variable = temp
    type = FunctionAux
    function = temp_function
  []
  [stress_zz_core]
    type = CylindricalRankTwoAux
    variable = stress_zz_core
    rank_two_tensor = stress
    index_j = 2
    index_i = 2
    center_point = '0 0 0'
    block = 1
  []
  [stress_rr_core]
    type = CylindricalRankTwoAux
    variable = stress_rr_core
    rank_two_tensor = stress
    index_j = 0
    index_i = 0
    center_point = '0 0 0'
    block = 1
  []
  [stress_tt_core]
    type = CylindricalRankTwoAux
    variable = stress_tt_core
    rank_two_tensor = stress
    index_j = 1
    index_i = 1
    center_point = '0 0 0'
    block = 1
  []
  [stress_zz_skin]
    type = CylindricalRankTwoAux
    variable = stress_zz_skin
    rank_two_tensor = stress
    index_j = 2
    index_i = 2
    center_point = '0 0 0'
    block = 2
  []
  [stress_rr_skin]
    type = CylindricalRankTwoAux
    variable = stress_rr_skin
    rank_two_tensor = stress
    index_j = 0
    index_i = 0
    center_point = '0 0 0'
    block = 2
  []
  [stress_tt_skin]
    type = CylindricalRankTwoAux
    variable = stress_tt_skin
    rank_two_tensor = stress
    index_j = 1
    index_i = 1
    center_point = '0 0 0'
    block = 2
  []
  [eigenstrain_yy]
    type = RankTwoAux
    rank_two_tensor = eigenstrain
    variable = eigenstrain_yy
    index_i = 1
    index_j = 1
    execute_on = 'initial timestep_end'
  []
  [eigenstrain_xx]
    type = RankTwoAux
    rank_two_tensor = eigenstrain
    variable = eigenstrain_xx
    index_i = 0
    index_j = 0
    execute_on = 'initial timestep_end'
  []
  [eigenstrain_zz]
    type = RankTwoAux
    rank_two_tensor = eigenstrain
    variable = eigenstrain_zz
    index_i = 2
    index_j = 2
    execute_on = 'initial timestep_end'
  []
[]
[Postprocessors]
  [temperature]
    type = ElementAverageValue
    variable = temp
  []
  [sigma_xx]
    type = ElementAverageValue
    variable = stress_xx
    execute_on = 'initial timestep_end'
  []
  [sigma_yy]
    type = ElementAverageValue
    variable = stress_yy
    execute_on = 'initial timestep_end'
  []
  [sigma_zz]
    type = ElementAverageValue
    variable = stress_zz
    execute_on = 'initial timestep_end'
  []
  [eigenstrain_xx]
    type = ElementAverageValue
    variable = eigenstrain_xx
    execute_on = 'initial timestep_end'
  []
  [eigenstrain_yy]
    type = ElementAverageValue
    variable = eigenstrain_yy
    execute_on = 'initial timestep_end'
  []
  [eigenstrain_zz]
    type = ElementAverageValue
    variable = eigenstrain_zz
    execute_on = 'initial timestep_end'
  []
[]
[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]
[Executioner]
  type = Transient
  solve_type = Newton
  start_time = 0.0
  dt = 0.1
  dtmax = 1.0
  dtmin = 0.001
  nl_abs_tol = 4e-07
  end_time = 1.0
  l_tol = 1e-5
[]
[Outputs]
  print_linear_residuals = true
  perf_graph = true
  csv = true
  [out]
    type = Exodus
    elemental_as_nodal = true
  []
  [console]
    type = Console
    max_rows = 100
  []
[]

Now the problem is if I change the Z-constraint to a single point: basically fix the origin [0.0, 0.0, 0.0] to Uz=0, the simulation no longer converge, which I believe is due to high stress generation at the constrained node.

However, if I impose no Uz constraints and only impose constraints in Ux and Uy, the simulation still runs giving a symmetric solution in stress which looks like below:

Now, my question is for the latter simulation where I did not impose any constraint along Z-axis how does MOOSE handle rigid body translation ?
And is there some numerical strategy which allows to solve for point constraint as boundary condition ?

Thanks,
Aritra

[GiudGiud]

Now, my question is for the latter simulation where I did not impose any constraint along Z-axis how does MOOSE handle rigid body translation ?

MOOSE does not handle it in any special way unless you add parameters to PETSc to remove the null space. (not the recommended solution, the recommended solution is to have the right constraints on the system to prevent RBM)
Without that, you are getting whatever the nonlinear solver returned, which may include rigid body motion, and which likely depends on the preconditioning method used and the number of solve iterations.

[GiudGiud]

is the linear solve diverging? You did not set petsc options so you are using the default preconditioning (ilu)
or is it the nonlinear solve?

[chakra34]

The linear residual seem to get stuck at a fixed value, and yeah I did not set any petsc options.

[chakra34]

[GiudGiud]

try setting

  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'

[chakra34]

damn! it solved!!

what's the magic recipe :D ?