PhaseField; energy split model.

Previously only the Volumetric-Deviatoric model was available. Now an Isotropic model is added.

MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h

@@ -32,10 +32,9 @@ std::tuple<MathLib::KelvinVector::KelvinVectorType<
                DisplacementDim> /* C_tensile */,
            MathLib::KelvinVector::KelvinMatrixType<
                DisplacementDim> /* C_compressive */,
-           double /* strain_energy_tensile */,
-           double /* elastic_energy */
+           double /* strain_energy_tensile */, double /* elastic_energy */
            >
-calculateDegradedStress(
+calculateVolDevDegradedStress(
     double const degradation,
     double const bulk_modulus,
     double const mu,
@@ -89,6 +88,56 @@ calculateDegradedStress(
     return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
                            strain_energy_tensile, elastic_energy);
 }
+
+template <int DisplacementDim>
+std::tuple<MathLib::KelvinVector::KelvinVectorType<
+               DisplacementDim> /* sigma_real */,
+           MathLib::KelvinVector::KelvinVectorType<
+               DisplacementDim> /* sigma_tensile */,
+           MathLib::KelvinVector::KelvinMatrixType<
+               DisplacementDim> /* C_tensile */,
+           MathLib::KelvinVector::KelvinMatrixType<
+               DisplacementDim> /* C_compressive */,
+           double /* strain_energy_tensile */, double /* elastic_energy */
+           >
+calculateIsotropicDegradedStress(
+    double const degradation,
+    double const bulk_modulus,
+    double const mu,
+    MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps)
+{
+    static int const KelvinVectorSize =
+        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+    using KelvinVector =
+        MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+    using KelvinMatrix =
+        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
+    using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+    // calculation of deviatoric parts
+    auto const& P_dev = Invariants::deviatoric_projection;
+    KelvinVector const epsd_curr = P_dev * eps;
+
+    // Hydrostatic part for the stress and the tangent.
+    double const eps_curr_trace = Invariants::trace(eps);
+
+    KelvinMatrix C_tensile = KelvinMatrix::Zero();
+    KelvinMatrix C_compressive = KelvinMatrix::Zero();
+
+    double const strain_energy_tensile =
+        bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
+        mu * epsd_curr.transpose() * epsd_curr;
+    double const elastic_energy = degradation * strain_energy_tensile;
+    KelvinVector const sigma_tensile =
+        bulk_modulus * eps_curr_trace * Invariants::identity2 +
+        2 * mu * epsd_curr;
+    C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
+    C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
+    KelvinVector const sigma_real = degradation * sigma_tensile;
+
+    return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
+                           strain_energy_tensile, elastic_energy);
+}
+
 }  // namespace Phasefield
 }  // namespace Solids
 }  // namespace MaterialLib

ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp

@@ -195,13 +195,33 @@ std::unique_ptr<Process> createPhaseFieldProcess(
             phasefield_model_string.c_str());
     }();
 
+    auto const energy_split_model = [&]
+    {
+        auto const energy_split_model_string =
+            //! \ogs_file_param{prj__processes__process__PHASE_FIELD__energy_split_model}
+            config.getConfigParameter<std::string>("energy_split_model");
+
+        if (energy_split_model_string == "Isotropic")
+        {
+            return EnergySplitModel::Isotropic;
+        }
+        else if (energy_split_model_string == "VolumetricDeviatoric")
+        {
+            return EnergySplitModel::VolDev;
+        }
+        OGS_FATAL(
+            "energy_split_model must be 'Isotropic' or 'VolumetricDeviatoric' "
+            "but '{:s}' was given",
+            energy_split_model_string);
+    }();
+
     PhaseFieldProcessData<DisplacementDim> process_data{
         materialIDs(mesh),   std::move(solid_constitutive_relations),
         residual_stiffness,  crack_resistance,
         crack_length_scale,  solid_density,
         specific_body_force, hydro_crack,
         crack_pressure,      irreversible_threshold,
-        phasefield_model};
+        phasefield_model,    energy_split_model};
 
     SecondaryVariableCollection secondary_variables;
 

ProcessLib/PhaseField/PhaseFieldFEM-impl.h

@@ -95,8 +95,9 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         double const k = _process_data.residual_stiffness(t, x_position)[0];
         double const d_ip = N.dot(d);
         double const degradation = d_ip * d_ip * (1 - k) + k;
-        _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
-                                                degradation);
+        _ip_data[ip].updateConstitutiveRelation(
+            t, x_position, dt, u, degradation,
+            _process_data.energy_split_model);
 
         auto& sigma = _ip_data[ip].sigma;
         auto const& C_tensile = _ip_data[ip].C_tensile;
@@ -184,8 +185,9 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
 
         auto& eps = _ip_data[ip].eps;
         eps.noalias() = B * u;
-        _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
-                                                degradation);
+        _ip_data[ip].updateConstitutiveRelation(
+            t, x_position, dt, u, degradation,
+            _process_data.energy_split_model);
 
         auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
 

ProcessLib/PhaseField/PhaseFieldFEM.h

@@ -74,7 +74,8 @@ struct IntegrationPointData final
                                     ParameterLib::SpatialPosition const& x,
                                     double const /*dt*/,
                                     DisplacementVectorType const& /*u*/,
-                                    double const degradation)
+                                    double const degradation,
+                                    EnergySplitModel const energy_split_model)
     {
         auto linear_elastic_mp =
             static_cast<MaterialLib::Solids::LinearElasticIsotropic<
@@ -84,22 +85,35 @@ struct IntegrationPointData final
         auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
         auto const mu = linear_elastic_mp.mu(t, x);
 
-        std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
-                 strain_energy_tensile, elastic_energy) =
-            MaterialLib::Solids::Phasefield::calculateDegradedStress<
-                DisplacementDim>(degradation, bulk_modulus, mu, eps);
+        switch (energy_split_model)
+        {
+            case EnergySplitModel::Isotropic:
+            {
+                std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
+                         strain_energy_tensile, elastic_energy) =
+                    MaterialLib::Solids::Phasefield::
+                        calculateIsotropicDegradedStress<DisplacementDim>(
+                            degradation, bulk_modulus, mu, eps);
+                break;
+            }
+            case EnergySplitModel::VolDev:
+            {
+                std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
+                         strain_energy_tensile, elastic_energy) =
+                    MaterialLib::Solids::Phasefield::
+                        calculateVolDevDegradedStress<DisplacementDim>(
+                            degradation, bulk_modulus, mu, eps);
+                break;
+            }
+        }
 
         history_variable =
             std::max(history_variable_prev, strain_energy_tensile);
     }
-    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
-
-    static constexpr int kelvin_vector_size =
-        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
 };
 
-/// Used for the extrapolation of the integration point values. It is ordered
-/// (and stored) by integration points.
+/// Used for the extrapolation of the integration point values. It is
+/// ordered (and stored) by integration points.
 template <typename ShapeMatrixType>
 struct SecondaryData
 {
@@ -242,7 +256,8 @@ public:
     }
 
     void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
-                              double const /*t*/, double const /*dt*/) override
+                              double const /*t*/,
+                              double const /*dt*/) override
     {
         unsigned const n_integration_points =
             _integration_method.getNumberOfPoints();

ProcessLib/PhaseField/PhaseFieldProcessData.h

@@ -38,6 +38,12 @@ enum class PhaseFieldModel
     AT2
 };
 
+enum class EnergySplitModel
+{
+    Isotropic,
+    VolDev
+};
+
 template <int DisplacementDim>
 struct PhaseFieldProcessData
 {
@@ -55,6 +61,7 @@ struct PhaseFieldProcessData
     bool crack_pressure = false;
     double irreversible_threshold;
     PhaseFieldModel phasefield_model;
+    EnergySplitModel energy_split_model;
 
     double const unity_pressure = 1.0;
     double pressure = 0.0;

ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h

@@ -92,7 +92,7 @@ struct IntegrationPointData final
 
         std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
                  strain_energy_tensile, elastic_energy) =
-            MaterialLib::Solids::Phasefield::calculateDegradedStress<
+            MaterialLib::Solids::Phasefield::calculateVolDevDegradedStress<
                 DisplacementDim>(degradation, bulk_modulus, mu, eps_m);
     }
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;