Variable coupled ElasticityTensor and Eigenstrain calculation

[keassyguang]

Hello everyone,

I am trying to develop a model to calculate a thermo-mechanical problem. I have successfully compiled the code, but the simulation has aborted with a "Segmentation fault (signal 11)" message. I have also tried to debug the code using GDB (example 21), but I haven't got any useful information.

The total strain is expressed as: ϵtotal = ϵelastic + ϵthermal + ϵchemical.

For thermal strain: ϵthermal = CTE(T-Tref), where CTE = cte_curedp + cte_uncured(1-p), where cte_cured and cte_uncured are two constants, and p is an ADMaterialProperty variable.

For chemical strain: ϵchemical = ϵchemical_max*p, where ϵchemical_max is a constant.

The Young's modulus varies with the ADMaterialProperty variable p, the relationship can be expressed as: E = Ecexp(bp), where Ec is the base Young's modulus and b is the exponential factor.

I have developed two classes: CureDependentIsotropicElasticity (based on #22587) and ChemicalThermalEigenStrain. Can anyone help me check my code to see if anything is wrong? Many thanks.

CureDependentIsotropicElasticity.h

#pragma once

#include "ComputeElasticityTensorBase.h"

/**
 * Compute an elasticity tensor that depends on the degree of cure
 * E = Ec * exp(b * p), where p is the degree of cure
 */

class CureDependentIsotropicElasticity : public ADComputeElasticityTensorBase
{
public:
  static InputParameters validParams();
  CureDependentIsotropicElasticity(const InputParameters & parameters);

protected:
  virtual void computeQpElasticityTensor() override;

  /// Base Young's modulus (Ec)
  const Real _Ec;
  
  /// Exponential factor for cure-dependent modulus (b)
  const Real _b;
  
  /// Poisson's ratio
  ADReal _youngs_modulus;
  ADReal _poissons_ratio;
  
  /// Degree of cure
  const ADMaterialProperty<Real> & _p;

  using ADComputeElasticityTensorBase::issueGuarantee;
  using ADComputeElasticityTensorBase::_elasticity_tensor;
  using ADComputeElasticityTensorBase::_qp;
};

CureDependentIsotropicElasticity.C

#include "CureDependentIsotropicElasticity.h"

registerMooseObject("MultiLayerDLPPolymerizationThermalMechanicalApp", CureDependentIsotropicElasticity);

InputParameters
CureDependentIsotropicElasticity::validParams()
{
  InputParameters params = ADComputeElasticityTensorBase::validParams();
  params.addParam<Real>("Ec", 0.08833e6, "Base Young's modulus (Ec)");
  params.addParam<Real>("b", 10.45, "Exponential factor (b)");
  params.addRequiredParam<Real>("poissons_ratio", "Poisson's ratio");
  params.addRequiredParam<MaterialPropertyName>("p", "Degree of cure property name");
  return params;
}

CureDependentIsotropicElasticity::CureDependentIsotropicElasticity(
    const InputParameters & parameters)
  : ADComputeElasticityTensorBase(parameters),
    _Ec(getParam<Real>("Ec")),
    _b(getParam<Real>("b")),
    _poissons_ratio(getParam<Real>("poissons_ratio")),
    _p(getADMaterialProperty<Real>("p"))
{
  // all tensors created by this class are always isotropic
  issueGuarantee(_elasticity_tensor_name, Guarantee::ISOTROPIC);
}

void
CureDependentIsotropicElasticity::computeQpElasticityTensor()
{
  _youngs_modulus = _Ec * std::exp(_b * _p[_qp]);

  // Assign elasticity tensor at a given quad point
  _elasticity_tensor[_qp].fillSymmetricIsotropicEandNu(_youngs_modulus, _poissons_ratio);
}

ChemicalThermalEigenStrain.h

#pragma once

#include "ADComputeEigenstrainBase.h"

/**
 * Compute eigenstrains due to both thermal expansion and chemical shrinkage
 */

class ChemicalThermalEigenStrain : public ADComputeEigenstrainBase
{
public:
  static InputParameters validParams();
  ChemicalThermalEigenStrain(const InputParameters & parameters);

protected:
  virtual void computeQpEigenstrain() override;

  /// Temperature variable
  const ADVariableValue & _temperature;
  
  /// Reference temperature
  const Real _reference_temperature;
  
  /// Degree of cure
  const ADMaterialProperty<Real> & _p;
  
  /// Maximum chemical shrinkage strain
  const Real _chemical_strain_max;
  
  /// Thermal expansion coefficient for cured material
  const Real _cte_cured;
  
  /// Thermal expansion coefficient for uncured material
  const Real _cte_uncured;
};

ChemicalThermalEigenStrain.C

#include "ChemicalThermalEigenStrain.h"

registerMooseObject("MultiLayerDLPPolymerizationThermalMechanicalApp", ChemicalThermalEigenStrain);

InputParameters
ChemicalThermalEigenStrain::validParams()
{
  InputParameters params = ADComputeEigenstrainBase::validParams();
  params.addRequiredCoupledVar("temperature", "Temperature variable");
  params.addRequiredParam<Real>("reference_temperature", "Reference temperature");
  params.addRequiredParam<MaterialPropertyName>("p", "Degree of cure property name");
  params.addParam<Real>("chemical_strain_max", -0.08, "Maximum chemical shrinkage strain");
  params.addParam<Real>("cte_cured", 50.0e-6, "Thermal expansion coefficient for cured material");
  params.addParam<Real>("cte_uncured", 156.0e-6, "Thermal expansion coefficient for uncured material");
  return params;
}

ChemicalThermalEigenStrain::ChemicalThermalEigenStrain(const InputParameters & parameters)
  : ADComputeEigenstrainBase(parameters),
    _temperature(adCoupledValue("temperature")),
    _reference_temperature(getParam<Real>("reference_temperature")),
    _p(getADMaterialProperty<Real>("p")),
    _chemical_strain_max(getParam<Real>("chemical_strain_max")),
    _cte_cured(getParam<Real>("cte_cured")),
    _cte_uncured(getParam<Real>("cte_uncured"))
{
}

void
ChemicalThermalEigenStrain::computeQpEigenstrain()
{
  ADReal cte = _p[_qp] * _cte_cured + (1.0 - _p[_qp]) * _cte_uncured;
  ADReal thermal_strain = cte * (_temperature[_qp] - _reference_temperature);
  ADReal chemical_strain = _chemical_strain_max * _p[_qp];
  
  _eigenstrain[_qp].zero();
  _eigenstrain[_qp].addIa(thermal_strain + chemical_strain);
}

Part of the input file is as follow:

[Variables]
  [disp_x]
    block = 1
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    block = 1
    order = FIRST
    family = LAGRANGE
  []
[]

[Kernels]
  # Mechanics kernels
  [stress_x]
    type = ADStressDivergenceTensors
    block = 1
    displacements = 'disp_x disp_y'
    variable = disp_x
    component = 0
    use_displaced_mesh = true
  []
  [stress_y]
    type = ADStressDivergenceTensors
    block = 1
    displacements = 'disp_x disp_y'
    variable = disp_y
    component = 1
    use_displaced_mesh = true
  []

[Materials]
  [elasticity]
    type = CureDependentIsotropicElasticity
    block = 1
    p = p
    poissons_ratio = 0.4
  []

  [strain]
    type = ADComputeFiniteStrain
    block = 1
    displacements = 'disp_x disp_y'
  []

  [eigenstrain]
    type = ChemicalThermalEigenStrain
    block = 1
    temperature = Temperature
    reference_temperature = 300.0
    p = p
    eigenstrain_name = chemical_thermal
  []

  [stress]
    type = ADComputeFiniteStrainElasticStress
    block = 1
  []
[]

[BCs]
  # BCs for mechanics
  # Fixed bottom boundary
  [bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'bottom'
    value = 0.0
  []
  [bottom_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
[]

[GiudGiud]

hello

but the simulation has aborted with a "Segmentation fault (signal 11)" message. I have also tried to debug the code using GDB (example 21), but I haven't got any useful information.

Can you report the backtrace please? This should be possible on a segfault

[keassyguang]

I still haven't solved the adaptive segmentation error problem. As I mentioned, the adaptivity works well with my thermal model, but it fails when the mechanics are introduced. I have checked the mesh, the mesh at the 'bottom' boundary starts to change at step 3. Do you have any suggestions to solve this problem?

Step 1: time = 0.1

Step 2: time = 0.2

Step 3: time = 0.3

[GiudGiud]

What do you mean by it fails?

Does it stop converging?

The Physics syntax is recommended over kernels as it's shorter and prevents some errors

[GiudGiud]

Oh ok it seg faults.

Let s a make a minimal non working example of this and I'll investigate it

[keassyguang]

What do you mean by it fails?

The simulation aborted with an error message.

Execution Information:
  Executioner:             Transient
  TimeStepper:             ConstantDT
  TimeIntegrator(s):       ImplicitEuler
  Solver Mode:             Preconditioned JFNK
  PETSc Preconditioner:    lu 


Time Step 0, time = 0

Time Step 1, time = 0.1, dt = 0.1

Performing automatic scaling calculation

 Solve Converged!

Time Step 2, time = 0.2, dt = 0.1
 Solve Converged!
Mesh unchanged, skipping remaining steps...

Time Step 3, time = 0.3, dt = 0.1
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Arguments are incompatible
[0]PETSC ERROR: Incompatible vector local lengths parameter # 1 local size 14440 != parameter # 2 local size 12546
[0]PETSC ERROR: WARNING! There are unused option(s) set! Could be the program crashed before usage or a spelling mistake, etc!
[0]PETSC ERROR:   Option left: name:-i value: test_multi_layer_polymerization_thermal_mechanical_with_AMR.i source: command line
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: unknown  GIT Date: unknown
[0]PETSC ERROR: ../multi_layer_dlp_polymerization_thermal_mechanical-opt with 1 MPI process(es) and PETSC_ARCH  on IHP-SPD-C0QH6A1 by liug2 Fri Jan 17 19:52:12 2025

The Physics syntax is recommended over kernels as it's shorter and prevents some errors

Thank you very much.

[keassyguang]

Let s a make a minimal non working example of this and I'll investigate it

Thank you very much, but it is not easy to find a simple example that does not work. I have tried to modify some of the examples provided by MOOSE, but all of them can be run successfully. Since we know that the problems might come from the adaptivity, I decide to try different parameters in the adaptivity section. There are several findings:

1. I tried different max_h_level parameters from 4 to 6, this problem seems to work well with max_h_level = 5, even with '-pc_type = hypre';
2. I tried to add another marker with a variable of von_mises, by carefully setting the threshold, the simulation can run successfully. With the same parameter, if the block was restricted to block 1, the simulation failed;
3. I also tried using a larger nl_abs_tol (1e-10 --> 1e-8) but it seems have no effects on the final solution.

As the simulation develops three different physics, the order of magnitude of the variables could be completely different, it seems that choosing a suitable indicator and marker to use the AMR certainly have significant effects on the convergence. I have seen some examples from malamute using both temperature and displacement variables as indicators. Do you have any recommendations on how to set an appropriate indicator and marker?

Beside, I also raise another discussion regarding the use of displaced mesh How to ensure the displaced mesh is used #29699. Can you offer me some solutions? Thanks a lot.