Merge pull request idaholab#311 from albiedahal/vectordamage

Vector Damage closes idaholab#307

include/materials/SteelCreepDamageOh.h

@@ -33,11 +33,14 @@ class SteelCreepDamageOhTempl : public ScalarDamageBaseTempl<is_ad>, public Guar
 protected:
   virtual void updateQpDamageIndex() override;
 
-  std::string _creep_strain_name;
+  std::vector<std::string> _creep_strain_names;
 
-  ///@{ Creep strain tensors
-  const GenericMaterialProperty<RankTwoTensor, is_ad> & _creep_strain;
-  const MaterialProperty<RankTwoTensor> & _creep_strain_old;
+  /// Creep strain tensors for each creep model
+  std::vector<const GenericMaterialProperty<RankTwoTensor, is_ad> *> _creep_model;
+
+  ///@{ Combined creep strain from all creep models
+  GenericMaterialProperty<RankTwoTensor, is_ad> & _combined_creep_strain;
+  const MaterialProperty<RankTwoTensor> & _combined_creep_strain_old;
   ///@}
 
   /// Uniaxial creep fracture strain (used to compute multiaxial creep fracture strain)

src/materials/SteelCreepDamageOh.C

@@ -44,21 +44,21 @@ SteelCreepDamageOhTempl<is_ad>::validParams()
       "reduction_damage_threshold <= 1.0 & reduction_damage_threshold >= 0.0",
       "Starting value of damage that will trigger linear reduction of "
       "quadrature point's load-carrying capacity.");
-  params.addParam<std::string>("creep_strain_name",
-                               "creep_strain",
-                               "Name of the creep strain material driving damage failure.");
+  params.addRequiredParam<std::vector<std::string>>(
+      "creep_strain_names", "Names of the creep strains driving the damage.");
   return params;
 }
 
 template <bool is_ad>
 SteelCreepDamageOhTempl<is_ad>::SteelCreepDamageOhTempl(const InputParameters & parameters)
   : ScalarDamageBaseTempl<is_ad>(parameters),
     GuaranteeConsumer(this),
-    _creep_strain_name(this->template getParam<std::string>("creep_strain_name")),
-    _creep_strain(this->template getGenericMaterialProperty<RankTwoTensor, is_ad>(
-        _base_name + _creep_strain_name)),
-    _creep_strain_old(
-        this->template getMaterialPropertyOld<RankTwoTensor>(_base_name + _creep_strain_name)),
+    _creep_strain_names(this->template getParam<std::vector<std::string>>("creep_strain_names")),
+    _creep_model(_creep_strain_names.size()),
+    _combined_creep_strain(this->template declareGenericProperty<RankTwoTensor, is_ad>(
+        _base_name + "combined_creep_strain")),
+    _combined_creep_strain_old(
+        this->template getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_creep_strain")),
     _epsilon_f(this->template getParam<Real>("epsilon_f")),
     _reduction_factor(this->template getParam<Real>("reduction_factor")),
     _reduction_damage_threshold(this->template getParam<Real>("reduction_damage_threshold")),
@@ -67,6 +67,11 @@ SteelCreepDamageOhTempl<is_ad>::SteelCreepDamageOhTempl(const InputParameters &
     _omega(this->template declareGenericProperty<Real, is_ad>(_base_name + "omega")),
     _omega_old(this->template getMaterialPropertyOld<Real>(_base_name + "omega"))
 {
+  for (unsigned int i = 0; i < _creep_strain_names.size(); ++i)
+  {
+    _creep_model[i] = &this->template getGenericMaterialProperty<RankTwoTensor, is_ad>(
+        _base_name + _creep_strain_names[i]);
+  }
   if (MooseUtils::absoluteFuzzyEqual(_creep_law_exponent, -0.5, TOLERANCE))
     this->template paramError(
         "creep_law_exponent",
@@ -80,6 +85,7 @@ SteelCreepDamageOhTempl<is_ad>::initQpStatefulProperties()
 {
   ScalarDamageBaseTempl<is_ad>::initQpStatefulProperties();
   _omega[_qp] = 0.0;
+  _combined_creep_strain[_qp].zero();
 }
 
 template <bool is_ad>
@@ -117,7 +123,12 @@ SteelCreepDamageOhTempl<is_ad>::updateQpDamageIndex()
     return;
   }
 
-  GenericRankTwoTensor<is_ad> creep_increment = _creep_strain[_qp] - _creep_strain_old[_qp];
+  _combined_creep_strain[_qp].zero();
+  for (unsigned int i = 0; i < _creep_strain_names.size(); ++i)
+    _combined_creep_strain[_qp] += (*_creep_model[i])[_qp];
+
+  GenericRankTwoTensor<is_ad> creep_increment =
+      _combined_creep_strain[_qp] - _combined_creep_strain_old[_qp];
 
   // Avoid derivative's divide by zero error
   Real epsilon_ad = 1.0e-14;

test/tests/steel_creep_damage_oh/ad_steel_creep_damage.i

@@ -120,7 +120,7 @@
   [damage]
     type = ADSteelCreepDamageOh
     epsilon_f = 0.01
-    creep_strain_name = creep_strain
+    creep_strain_names = creep_strain
     reduction_factor = 1.0e3
     use_old_damage = true
     creep_law_exponent = 10.0

test/tests/steel_creep_damage_oh/ad_steel_creep_damage_2d.i

@@ -113,7 +113,7 @@
   [damage]
     type = ADSteelCreepDamageOh
     epsilon_f = 0.01
-    creep_strain_name = creep_strain
+    creep_strain_names = creep_strain
     reduction_factor = 1.0e3
     use_old_damage = true
     creep_law_exponent = 10.0

test/tests/steel_creep_damage_oh/gold/multi_power_law_creep_out.csv

@@ -0,0 +1 @@
+power_law_creep_out.csv

test/tests/steel_creep_damage_oh/gold/power_law_creep_out.csv

@@ -0,0 +1,7 @@
+time,combined_creep_strain_yy,damage_index,dt,max_disp_x,max_disp_y,omega
+0,0,0,0,0,0,0
+0.025,1.3591379713935e-07,0,0.025,-6.3775703573435e-08,2.3591229183113e-06,5.1088169606264e-06
+0.05,1.1224595523922e-06,0,0.025,-1.5009847904727e-07,4.7182458365616e-06,4.2517326705618e-05
+0.075,2.9728643440701e-06,0,0.025,-2.6193875410219e-07,7.0773687552007e-06,0.00011744316468622
+0.1,5.2481408338636e-06,0,0.025,-3.8859149372669e-07,9.4364916740147e-06,0.00021475605982709
+0.125,7.6822935439709e-06,0,0.025,-5.2230942992542e-07,1.1795614592695e-05,0.00032310458677235

test/tests/steel_creep_damage_oh/power_law_creep.i

@@ -0,0 +1,187 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  nx = 2
+  ny = 2
+  second_order = true
+[]
+
+[GlobalParams]
+  displacements = 'disp_x disp_y'
+  volumetric_locking_correction = false
+[]
+
+[AuxVariables]
+  [local_damage]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [omega]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [combined_creep_strain_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+[]
+
+[AuxKernels]
+  [local_damage]
+    type = MaterialRealAux
+    variable = local_damage
+    property = local_damage
+    execute_on = TIMESTEP_END
+  []
+  [omega]
+    type = MaterialRealAux
+    variable = omega
+    property = omega
+  []
+  [combined_creep_strain_yy]
+    type = RankTwoAux
+    variable = combined_creep_strain_yy
+    rank_two_tensor = combined_creep_strain
+    index_j = 1
+    index_i = 1
+    execute_on = timestep_end
+  []
+[]
+
+[Functions]
+  [pull]
+    type = PiecewiseLinear
+    x = '0 10'
+    y = '0 1e-3'
+  []
+[]
+
+[Modules/TensorMechanics/Master]
+  [all]
+    strain = FINITE
+    add_variables = true
+    use_finite_deform_jacobian = true
+    generate_output = 'hydrostatic_stress'
+  []
+[]
+
+[Materials]
+  [elasticity_tensor]
+    type = ComputeIsotropicElasticityTensor
+    youngs_modulus = 1e10
+    poissons_ratio = 0.3
+  []
+
+  [elastic_strain]
+    type = ComputeMultipleInelasticStress
+    # tangent_operator = nonlinear
+    damage_model = 'steel_local'
+  []
+
+  [steel_local]
+    type = SteelCreepDamageOh
+    epsilon_f = 0.07
+    creep_law_exponent = 10
+    use_displaced_mesh = true
+    damage_index_name = local_damage
+    reduction_damage_threshold = 0.5
+    use_old_damage = true
+  []
+  [creep_ten]
+    type = PowerLawCreepStressUpdate
+    coefficient = 10e-24
+    n_exponent = 4
+    activation_energy = 0
+    base_name = creep_ten
+  []
+  [creep_ten2]
+    type = PowerLawCreepStressUpdate
+    coefficient = 10e-24
+    n_exponent = 4
+    activation_energy = 0
+    base_name = creep_ten2
+  []
+  # This model is activated with cli_args to compare the response with a single model with that of two models
+  [creep_twenty]
+    type = PowerLawCreepStressUpdate
+    coefficient = 20e-24
+    n_exponent = 4
+    activation_energy = 0
+    base_name = creep_twenty
+  []
+[]
+
+[BCs]
+  [no_disp_x]
+    type = DirichletBC
+    variable = disp_x
+    boundary = left
+    value = 0.0
+  []
+
+  [no_disp_y]
+    type = DirichletBC
+    variable = disp_y
+    boundary = bottom
+    value = 0.0
+  []
+
+  [pull_disp_y]
+    type = FunctionDirichletBC
+    variable = disp_y
+    boundary = top
+    function = pull
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = PJFNK
+
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+  line_search = 'none'
+  nl_rel_tol = 1e-10
+  nl_abs_tol = 1e-50
+
+  num_steps = 5
+  dt = 0.025
+[]
+
+[Postprocessors]
+  [max_disp_x]
+    type = ElementExtremeValue
+    variable = disp_x
+  []
+  [max_disp_y]
+    type = ElementExtremeValue
+    variable = disp_y
+  []
+  [dt]
+    type = TimestepSize
+  []
+  [damage_index]
+    type = ElementExtremeValue
+    variable = local_damage
+  []
+  [omega]
+    type = ElementAverageValue
+    variable = omega
+  []
+  [combined_creep_strain_yy]
+    type = ElementAverageValue
+    variable = combined_creep_strain_yy
+  []
+[]
+
+[Outputs]
+  csv = true
+  perf_graph = true
+[]

test/tests/steel_creep_damage_oh/steel_creep_damage.i

@@ -119,7 +119,7 @@
   [damage]
     type = SteelCreepDamageOh
     epsilon_f = 0.01
-    creep_strain_name = creep_strain
+    creep_strain_names = creep_strain
     reduction_factor = 1.0e3
     use_old_damage = true
     creep_law_exponent = 10.0

test/tests/steel_creep_damage_oh/steel_creep_damage_2d.i

@@ -112,7 +112,7 @@
   [damage]
     type = SteelCreepDamageOh
     epsilon_f = 0.01
-    creep_strain_name = creep_strain
+    creep_strain_names = creep_strain
     reduction_factor = 1.0e3
     use_old_damage = true
     creep_law_exponent = 10.0