(M) Constitutive laws to consider dimensionless formulation

include/gCP/constitutive_laws.h

@@ -26,7 +26,9 @@ template <int dim>
 class ElasticStrain
 {
 public:
-  ElasticStrain(std::shared_ptr<CrystalsData<dim>>  crystals_data);
+  ElasticStrain(std::shared_ptr<CrystalsData<dim>>  crystals_data,
+  const double reference_length_value = 1.0,
+  const double reference_displacement_value = 1.0);
 
   const dealii::SymmetricTensor<2,dim> get_elastic_strain_tensor(
     const unsigned int                      crystal_id,
@@ -41,6 +43,8 @@ class ElasticStrain
 
 private:
   std::shared_ptr<const CrystalsData<dim>>    crystals_data;
+
+  double                                      dimensionless_number;
 };
 
 
@@ -57,10 +61,12 @@ template<int dim>
 class HookeLaw
 {
 public:
-  HookeLaw(const RunTimeParameters::HookeLawParameters  parameters);
+  HookeLaw(const RunTimeParameters::HookeLawParameters  parameters,
+           const double reference_stiffness_value = 1.0);
 
   HookeLaw(const std::shared_ptr<CrystalsData<dim>>     &crystals_data,
-           const RunTimeParameters::HookeLawParameters  parameters);
+           const RunTimeParameters::HookeLawParameters  parameters,
+           const double reference_stiffness_value = 1.0);
 
   void init();
 
@@ -106,6 +112,7 @@ class HookeLaw
 
   std::vector<dealii::SymmetricTensor<4,3>>   stiffness_tetrads_3d;
 
+  double                                      reference_stiffness_value;
 
   bool                                        flag_init_was_called;
 };
@@ -310,10 +317,10 @@ class VectorialMicrostressLaw
 {
 public:
   VectorialMicrostressLaw(
-    const std::shared_ptr<CrystalsData<dim>>                      &crystals_data,
+    const std::shared_ptr<CrystalsData<dim>> &crystals_data,
     const RunTimeParameters::VectorialMicrostressLawParameters parameters);
 
-  void init();
+  void init(const bool flag_dimensionless_formulation = false);
 
   dealii::Tensor<1,dim> get_vectorial_microstress(
     const unsigned int          crystal_id,
@@ -339,6 +346,8 @@ class VectorialMicrostressLaw
 
   const double                                defect_energy_index;
 
+  double                                      factor;
+
   std::vector<std::vector<dealii::SymmetricTensor<2,dim>>>
                                               slip_direction_dyads;
 

include/gCP/gradient_crystal_plasticity.h

@@ -111,6 +111,8 @@ class GradientCrystalPlasticitySolver
 
   std::shared_ptr<dealii::TensorFunction<1,dim>>    supply_term;
 
+  std::vector<double>                               dimensionless_numbers;
+
   std::shared_ptr<Kinematics::ElasticStrain<dim>>   elastic_strain;
 
   std::shared_ptr<ConstitutiveLaws::HookeLaw<dim>>  hooke_law;

include/gCP/quadrature_point_history.h

@@ -291,7 +291,8 @@ struct QuadraturePointHistory
    */
   void init(
     const RunTimeParameters::HardeningLaw &parameters,
-    const unsigned int                    n_slips);
+    const unsigned int                    n_slips,
+    const double reference_slip_resistance_value = 1.0);
 
   /*!
    * @brief Stores the values of @ref slip_resistances in @ref
@@ -346,6 +347,8 @@ struct QuadraturePointHistory
 
   double              hardening_parameter;
 
+  double              reference_slip_resistance_value;
+
   bool                flag_perfect_plasticity;
 
   bool                flag_init_was_called;

source/constitutive_laws.cc

@@ -16,9 +16,13 @@ namespace Kinematics
 
 template <int dim>
 ElasticStrain<dim>::ElasticStrain(
-  std::shared_ptr<CrystalsData<dim>>  crystals_data)
+  std::shared_ptr<CrystalsData<dim>>  crystals_data,
+  const double reference_length_value,
+  const double reference_displacement_value)
 :
-crystals_data(crystals_data)
+crystals_data(crystals_data),
+dimensionless_number(
+  reference_length_value / reference_displacement_value)
 {}
 
 
@@ -39,11 +43,11 @@ ElasticStrain<dim>::get_elastic_strain_tensor(
   dealii::SymmetricTensor<2,dim> elastic_strain_tensor_value(
                                   strain_tensor_value);
 
-  for (unsigned int slip_id = 0;
-       slip_id < crystals_data->get_n_slips();
-       ++slip_id)
+  for (unsigned int slip_id = 0; slip_id < crystals_data->get_n_slips();
+        ++slip_id)
   {
     elastic_strain_tensor_value -=
+      dimensionless_number *
       slip_values[slip_id][q_point] *
       crystals_data->get_symmetrized_schmid_tensor(crystal_id, slip_id);
   }
@@ -92,12 +96,14 @@ namespace ConstitutiveLaws
 
 template<int dim>
 HookeLaw<dim>::HookeLaw(
-  const RunTimeParameters::HookeLawParameters  parameters)
+  const RunTimeParameters::HookeLawParameters  parameters,
+  const double reference_stiffness_value)
 :
 crystallite(Crystallite::Monocrystalline),
 C1111(parameters.C1111),
 C1122(parameters.C1122),
 C1212(parameters.C1212),
+reference_stiffness_value(reference_stiffness_value),
 flag_init_was_called(false)
 {
   crystals_data = nullptr;
@@ -108,13 +114,15 @@ flag_init_was_called(false)
 template<int dim>
 HookeLaw<dim>::HookeLaw(
   const std::shared_ptr<CrystalsData<dim>>    &crystals_data,
-  const RunTimeParameters::HookeLawParameters parameters)
+  const RunTimeParameters::HookeLawParameters parameters,
+  const double reference_stiffness_value)
 :
 crystals_data(crystals_data),
 crystallite(Crystallite::Polycrystalline),
 C1111(parameters.C1111),
 C1122(parameters.C1122),
 C1212(parameters.C1212),
+reference_stiffness_value(reference_stiffness_value),
 flag_init_was_called(false)
 {}
 
@@ -128,11 +136,14 @@ void HookeLaw<dim>::init()
       for (unsigned int k = 0; k < dim; k++)
         for (unsigned int l = 0; l < dim; l++)
           if (i == j && j == k && k == l)
-            reference_stiffness_tetrad[i][j][k][l] = C1111;
+            reference_stiffness_tetrad[i][j][k][l] =
+              C1111 / reference_stiffness_value;
           else if (i == k && j == l)
-            reference_stiffness_tetrad[i][j][k][l] = C1212;
+            reference_stiffness_tetrad[i][j][k][l] =
+              C1212 / reference_stiffness_value;
           else if (i == j && k == l)
-            reference_stiffness_tetrad[i][j][k][l] = C1122;
+            reference_stiffness_tetrad[i][j][k][l] =
+              C1122 / reference_stiffness_value;
 
   if constexpr(dim == 3)
     reference_stiffness_tetrad_3d = reference_stiffness_tetrad;
@@ -143,11 +154,14 @@ void HookeLaw<dim>::init()
         for (unsigned int k = 0; k < 3; k++)
           for (unsigned int l = 0; l < 3; l++)
             if (i == j && j == k && k == l)
-              reference_stiffness_tetrad_3d[i][j][k][l] = C1111;
+              reference_stiffness_tetrad_3d[i][j][k][l] =
+                C1111 / reference_stiffness_value;
             else if (i == k && j == l)
-              reference_stiffness_tetrad_3d[i][j][k][l] = C1212;
+              reference_stiffness_tetrad_3d[i][j][k][l] =
+                C1212 / reference_stiffness_value;
             else if (i == j && k == l)
-              reference_stiffness_tetrad_3d[i][j][k][l] = C1122;
+              reference_stiffness_tetrad_3d[i][j][k][l] =
+                C1122 / reference_stiffness_value;
   }
   else
     Assert(false, dealii::ExcNotImplemented());
@@ -574,13 +588,19 @@ flag_init_was_called(false)
 
 
 template<int dim>
-void VectorialMicrostressLaw<dim>::init()
+void VectorialMicrostressLaw<dim>::init(
+  const bool flag_dimensionless_formulation)
 {
   AssertThrow(crystals_data->is_initialized(),
               dealii::ExcMessage("The underlying CrystalsData<dim>"
                                   " instance has not been "
                                   " initialized."));
 
+  factor = flag_dimensionless_formulation ?
+              1.0 :
+              initial_slip_resistance *
+                std::pow(energetic_length_scale, defect_energy_index);
+
   for (unsigned int crystal_id = 0;
         crystal_id < crystals_data->get_n_crystals();
         crystal_id++)
@@ -661,8 +681,7 @@ get_vectorial_microstress(
                                  "instance has not been initialized."));
 
   return (
-    initial_slip_resistance *
-    std::pow(energetic_length_scale, defect_energy_index) *
+    factor *
     (
       std::pow(
         std::abs(-crystals_data->get_slip_direction(crystal_id, slip_id) *
@@ -696,8 +715,7 @@ get_jacobian(
                                  "instance has not been initialized."));
 
   return (
-    initial_slip_resistance *
-    std::pow(energetic_length_scale, defect_energy_index) *
+    factor *
     (defect_energy_index - 1.0) *
     (
       std::pow(

source/gradient_crystal_plasticity/quadrature_point_history.cc

@@ -286,7 +286,8 @@ flag_init_was_called(false)
 template <int dim>
 void QuadraturePointHistory<dim>::init(
   const RunTimeParameters::HardeningLaw &parameters,
-  const unsigned int                    n_slips)
+  const unsigned int                    n_slips,
+  const double reference_slip_resistance_value)
 {
   this->n_slips             = n_slips;
 
@@ -300,10 +301,14 @@ void QuadraturePointHistory<dim>::init(
 
   slip_resistances          = std::vector<double>(
                                 n_slips,
-                                initial_slip_resistance);
+                                initial_slip_resistance /
+                                  reference_slip_resistance_value);
 
   tmp_slip_resistances      = slip_resistances;
 
+  this->reference_slip_resistance_value =
+    reference_slip_resistance_value;
+
   flag_init_was_called      = true;
 }
 
@@ -359,7 +364,8 @@ void QuadraturePointHistory<dim>::update_values(
           ++slip_id_beta)
     {
       slip_resistances[slip_id_alpha] +=
-        get_hardening_matrix_entry(slip_id_alpha == slip_id_beta) *
+        get_hardening_matrix_entry(slip_id_alpha == slip_id_beta) /
+        reference_slip_resistance_value *
         std::fabs(slips[slip_id_beta][q_point] -
                   old_slips[slip_id_beta][q_point]);
     }

source/gradient_crystal_plasticity/setup.cc

@@ -102,10 +102,32 @@ void GradientCrystalPlasticitySolver<dim>::init()
     jacobian.reinit(sparsity_pattern);
   }
 
+  // Compute reference parameters
+  dimensionless_numbers[0] =
+    parameters.reference_parameters.reference_length_value /
+      parameters.reference_parameters.reference_displacement_value;
+
+  dimensionless_numbers[1] =
+    parameters.constitutive_laws_parameters.hardening_law_parameters.
+      linear_hardening_modulus / parameters.reference_parameters.
+        reference_slip_resistance_value;
+
+  dimensionless_numbers[2] =
+    std::pow(parameters.constitutive_laws_parameters.
+      vectorial_microstress_law_parameters.energetic_length_scale /
+        parameters.reference_parameters.reference_length_value,
+          parameters.constitutive_laws_parameters.
+            vectorial_microstress_law_parameters.defect_energy_index);
+
+  dimensionless_numbers[3] =
+    parameters.reference_parameters.reference_stiffness_value /
+      parameters.reference_parameters.reference_slip_resistance_value /
+        dimensionless_numbers[0];
+
   // Initiate constitutive laws
   hooke_law->init();
 
-  vectorial_microstress_law->init();
+  vectorial_microstress_law->init(dimensionless_numbers[3] != 1.0);
 
   init_quadrature_point_history();
 
@@ -431,7 +453,9 @@ void GradientCrystalPlasticitySolver<dim>::init_quadrature_point_history()
             ++q_point)
         local_quadrature_point_history[q_point]->init(
           parameters.constitutive_laws_parameters.hardening_law_parameters,
-          crystals_data->get_n_slips());
+          crystals_data->get_n_slips(),
+          parameters.reference_parameters.
+            reference_slip_resistance_value);
 
       if (cell_is_at_grain_boundary(cell->active_cell_index()) &&
           fe_field->is_decohesion_allowed())