Phasefield (#1895)

* added phasefield process and test.

* apply clang-format.

* added missing integralMeasure for axisymmetric problem.

* fixed compilation errors after rebasing.

* added unconfined compression ctest.

* added missing reference.

* renamed specialFunction as calculateDegradedStress.

* changed order of stress and strain components.

* deleted unused file.

* changed IntegrationPointData.

* added local results to highlight the deviation.

* returned DeltaXCheck true.

* modified dimension in ProjectData.

* removed second order mesh.

* moved history_variable to pushBackState.

* added absolute error of displacement for beam3d test.

* removed eval() calls in calculation of stress components.

* added description for the calculateDegradedStress function.

* added types available.

* rephrased description.

* removed unused C variable.

* added model description.

* modified types of constitutive relation.

* removed CreateLinearElasticPhaseField.h.

* added a tag for phasefield_parameters.

* added missing phasefield_parameters after PHASE_FIELD.

* removed eval() calls in calculation of tensile stress component.

* moved process description to PhaseFieldProcess.h.

* added description to the undocumented function.

* created missing file in documentation.

* removed two duplicated files in phasefield process.

* removed unused file in documentation.

* removed the call to material_model.integrateStress.

* removed unused variables.

* fixed expression of length scale parameter in documentation.

Applications/ApplicationsLib/ProjectData.cpp

@@ -42,6 +42,7 @@
 #include "ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.h"
 #include "ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.h"
 #include "ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h"
+#include "ProcessLib/PhaseField/CreatePhaseFieldProcess.h"
 #include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
 #include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
 #include "ProcessLib/TES/CreateTESProcess.h"
@@ -399,6 +400,26 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
                     _process_variables, _parameters, integration_order,
                     process_config);
         }
+        else if (type == "PHASE_FIELD")
+        {
+            switch (_mesh_vec[0]->getDimension())
+            {
+                case 2:
+                    process =
+                        ProcessLib::PhaseField::createPhaseFieldProcess<
+                            2>(*_mesh_vec[0], std::move(jacobian_assembler),
+                               _process_variables, _parameters,
+                               integration_order, process_config);
+                    break;
+                case 3:
+                    process =
+                        ProcessLib::PhaseField::createPhaseFieldProcess<
+                            3>(*_mesh_vec[0], std::move(jacobian_assembler),
+                               _process_variables, _parameters,
+                               integration_order, process_config);
+                    break;
+            }
+        }
         else if (type == "SMALL_DEFORMATION")
         {
             switch (_mesh_vec[0]->getDimension())

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/c_PHASE_FIELD.md

+Phase-field modelling for brittle fracture. It is implemented monolithically.
\ No newline at end of file

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/constitutive_relation/i_constitutive_relation.md

+The constitutive relation for the brittle fracture.
\ No newline at end of file

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/constitutive_relation/t_type.md

+The type of constitutive relation available.
+1. Linear elastic
\ No newline at end of file

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/i_phasefield_parameters.md

+A tag for phase-field parameters.

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_crack_length_scale.md

+A length-scale parameter which controls the regularisation.
\ No newline at end of file

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_crack_resistance.md

+The critical Griffith-type fracture energy.
\ No newline at end of file

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_history_field.md

+A damage-driving history field which is associated with the maximum local tensile strain energy.
\ No newline at end of file

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_kinetic_coefficient.md

+A kinetic coefficient which is defined to address rate-dependent cases.
\ No newline at end of file

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_residual_stiffness.md

+A residual stiffness which is used in numerical calculations.

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/i_process_variables.md

+The process variables for displacement and phasefield.

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/t_displacement.md

+Process variable name of displacement.

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/t_phasefield.md

+Process variable name of phasefield.

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_solid_density.md

+Reference solid density.
\ No newline at end of file

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_specific_body_force.md

+Specific body forces that are used to apply gravitational forces. A vector of displacement dimension's length.

MaterialLib/SolidModels/LinearElasticIsotropic.h

@@ -17,7 +17,7 @@ namespace MaterialLib
 namespace Solids
 {
 template <int DisplacementDim>
-class LinearElasticIsotropic final : public MechanicsBase<DisplacementDim>
+class LinearElasticIsotropic : public MechanicsBase<DisplacementDim>
 {
 public:
     /// Variables specific to the material model
@@ -47,6 +47,13 @@ public:
                    (2 * (1 + _poissons_ratio(t, x)[0]));
         }
 
+        /// the bulk modulus.
+        double bulk_modulus(double const t, X const& x) const
+        {
+            return _youngs_modulus(t, x)[0] /
+                   (3 * (1 - 2 * _poissons_ratio(t, x)[0]));
+        }
+
     private:
         P const& _youngs_modulus;
         P const& _poissons_ratio;
@@ -101,7 +108,9 @@ public:
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
             material_state_variables) override;
 
-private:
+    MaterialProperties getMaterialProperties() {return _mp;}
+
+protected:
     MaterialProperties _mp;
 };
 

MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "LinearElasticIsotropicPhaseField.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+
+template class LinearElasticIsotropicPhaseField<2>;
+template class LinearElasticIsotropicPhaseField<3>;
+
+}   // namespace Solids
+}   // namespace MaterialLib

MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "PhaseFieldExtension.h"
+#include "LinearElasticIsotropic.h"
+#include "KelvinVector.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+class LinearElasticIsotropicPhaseField final
+    : public LinearElasticIsotropic<DisplacementDim>,
+      public PhaseFieldExtension<DisplacementDim>
+{
+public:
+    static int const KelvinVectorSize =
+        ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+    using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+    using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+
+    explicit LinearElasticIsotropicPhaseField(
+        typename LinearElasticIsotropic<
+            DisplacementDim>::MaterialProperties&& material_properties)
+        : LinearElasticIsotropic<DisplacementDim>(std::move(material_properties))
+    {
+    }
+
+    std::unique_ptr<
+        typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
+    createMaterialStateVariables() override
+    {
+        return LinearElasticIsotropic<
+            DisplacementDim>::createMaterialStateVariables();
+    }
+
+    boost::optional<std::tuple<KelvinVector,
+                               std::unique_ptr<typename MechanicsBase<
+                                   DisplacementDim>::MaterialStateVariables>,
+                               KelvinMatrix>>
+    integrateStress(
+        double const t,
+        ProcessLib::SpatialPosition const& x,
+        double const dt,
+        KelvinVector const& eps_prev,
+        KelvinVector const& eps,
+        KelvinVector const& sigma_prev,
+        typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
+            material_state_variables) override
+    {
+        return LinearElasticIsotropic<DisplacementDim>::
+                        integrateStress(t,
+                                        x,
+                                        dt,
+                                        eps_prev,
+                                        eps,
+                                        sigma_prev,
+                                        material_state_variables);
+    }
+
+    /** Decompose the stiffness into tensile and compressive part.
+     * Judging by the physical observations, compression perpendicular
+     * to a crack does not cause crack propagation. Thus,
+     * the phase-field parameter is only involved into the tensile part
+     * to degrade the elastic strain energy.
+     */
+    bool calculateDegradedStress(double const t,
+                         ProcessLib::SpatialPosition const& x,
+                         KelvinVector const& eps,
+                         double& strain_energy_tensile,
+                         KelvinVector& sigma_tensile,
+                         KelvinVector& sigma_compressive,
+                         KelvinMatrix& C_tensile,
+                         KelvinMatrix& C_compressive,
+                         KelvinVector& sigma_real,
+                         double const degradation) const override
+    {
+        using Invariants =
+            MaterialLib::SolidModels::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);
+
+        auto const& K =
+            LinearElasticIsotropic<DisplacementDim>::_mp.bulk_modulus(t, x);
+        auto const& mu =
+            LinearElasticIsotropic<DisplacementDim>::_mp.mu(t, x);
+
+        C_tensile = KelvinMatrix::Zero();
+        C_compressive = KelvinMatrix::Zero();
+
+        if (eps_curr_trace >= 0)
+        {
+            strain_energy_tensile =
+                 K / 2 * eps_curr_trace * eps_curr_trace +
+                 mu * epsd_curr.transpose() * epsd_curr;
+            sigma_tensile.noalias() =
+                 K * eps_curr_trace * Invariants::identity2 +
+                 2 * mu * epsd_curr;
+            sigma_compressive.noalias() = KelvinVector::Zero();
+            C_tensile.template topLeftCorner<3, 3>().setConstant(K);
+            C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
+        }
+        else
+        {
+            strain_energy_tensile = mu * epsd_curr.transpose() * epsd_curr;
+            sigma_tensile.noalias() = 2 * mu * epsd_curr;
+            sigma_compressive.noalias() =
+                 K * eps_curr_trace * Invariants::identity2;
+            C_tensile.noalias() = 2 * mu * P_dev * KelvinMatrix::Identity();
+            C_compressive.template topLeftCorner<3, 3>().setConstant(K);
+        }
+
+        sigma_real.noalias() = degradation * sigma_tensile + sigma_compressive;
+        return true;
+    }
+};
+
+extern template class LinearElasticIsotropicPhaseField<2>;
+extern template class LinearElasticIsotropicPhaseField<3>;
+
+}  // namespace Solids
+}  // namespace MaterialLib

MaterialLib/SolidModels/PhaseFieldExtension.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
+#include "MechanicsBase.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+
+template <int DisplacementDim>
+struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
+{
+    using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+    using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+    virtual bool calculateDegradedStress(double const t,
+                                 ProcessLib::SpatialPosition const& x,
+                                 KelvinVector const& eps,
+                                 double& strain_energy_tensile,
+                                 KelvinVector& sigma_tensile,
+                                 KelvinVector& sigma_compressive,
+                                 KelvinMatrix& C_tensile,
+                                 KelvinMatrix& C_compressive,
+                                 KelvinVector& sigma_real,
+                                 double const degradation) const = 0;
+
+    /// Dynamic size Kelvin vector and matrix wrapper for the polymorphic
+    /// constitutive relation compute function.
+    bool calculateDegradedStress(
+        double const t,
+        ProcessLib::SpatialPosition const& x,
+        Eigen::Matrix<double, Eigen::Dynamic, 1> const& eps,
+        double& strain_energy_tensile,
+        Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_tensile,
+        Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_compressive,
+        Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
+            C_tensile,
+        Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
+            C_compressive,
+        Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_real,
+        double const degradation) const
+    {
+        // TODO Avoid copies of data:
+        // Using MatrixBase<Derived> not possible because template functions
+        // cannot be virtual. Maybe there is a workaround for this.  Using
+        // Map<Matrix<double, ...>> makes the interface (for the material model
+        // implementation) unnecessary difficult.
+        KelvinVector const eps_{eps};
+        KelvinVector sigma_tensile_{sigma_tensile};
+        KelvinVector sigma_compressive_{sigma_compressive};
+        KelvinMatrix C_tensile_{C_tensile};
+        KelvinMatrix C_compressive_{C_compressive};
+        KelvinVector sigma_real_{sigma_real};
+
+        bool const result = calculateDegradedStress(t,
+                                            x,
+                                            eps_,
+                                            strain_energy_tensile,
+                                            sigma_tensile_,
+                                            sigma_compressive_,
+                                            C_tensile_,
+                                            C_compressive_,
+                                            sigma_real_,
+                                            degradation);
+
+        sigma_tensile = sigma_tensile_;
+        sigma_compressive = sigma_compressive_;
+        C_tensile = C_tensile_;
+        C_compressive = C_compressive_;
+        sigma_real = sigma_real_;
+        return result;
+    }
+
+    virtual ~PhaseFieldExtension() = default;
+};
+
+}  // namespace Solids
+}  // namespace MaterialLib

ProcessLib/CMakeLists.txt

@@ -17,6 +17,7 @@ APPEND_SOURCE_FILES(SOURCES LIE/SmallDeformation)
 APPEND_SOURCE_FILES(SOURCES LIE/SmallDeformation/LocalAssembler)
 APPEND_SOURCE_FILES(SOURCES LiquidFlow)
 APPEND_SOURCE_FILES(SOURCES Parameter)
+APPEND_SOURCE_FILES(SOURCES PhaseField)
 APPEND_SOURCE_FILES(SOURCES RichardsFlow)
 APPEND_SOURCE_FILES(SOURCES SmallDeformation)
 APPEND_SOURCE_FILES(SOURCES TES)
@@ -58,6 +59,7 @@ include(HydroMechanics/Tests.cmake)
 include(LIE/HydroMechanics/Tests.cmake)
 include(LIE/SmallDeformation/Tests.cmake)
 include(LiquidFlow/Tests.cmake)
+include(PhaseField/Tests.cmake)
 include(RichardsFlow/Tests.cmake)
 include(SmallDeformation/Tests.cmake)
 include(TES/Tests.cmake)