a degraded elasticity tensor in the weak form implemented

MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h

@@ -50,10 +50,7 @@ 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 */,
+           MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
            double /* strain_energy_tensile */, double /* elastic_energy */
            >
 calculateVolDevDegradedStress(
@@ -113,8 +110,11 @@ calculateVolDevDegradedStress(
         mu * epsd_curr.transpose() * epsd_curr;
     KelvinVector const sigma_real =
         degradation * sigma_tensile + sigma_compressive;
-    return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
-                           strain_energy_tensile, elastic_energy);
+
+    KelvinMatrix const D = degradation * C_tensile + C_compressive;
+
+    return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
+                           elastic_energy);
 }
 
 template <int DisplacementDim>
@@ -122,10 +122,7 @@ 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 */,
+           MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
            double /* strain_energy_tensile */, double /* elastic_energy */
            >
 calculateIsotropicDegradedStress(
@@ -162,8 +159,10 @@ calculateIsotropicDegradedStress(
     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);
+    KelvinMatrix const D = degradation * C_tensile + C_compressive;
+
+    return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
+                           elastic_energy);
 }
 
 template <int DisplacementDim>
@@ -171,10 +170,7 @@ 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 */,
+           MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
            double /* strain_energy_tensile */, double /* elastic_energy */
            >
 calculateIsotropicDegradedStressWithRankineEnergy(
@@ -227,8 +223,9 @@ calculateIsotropicDegradedStressWithRankineEnergy(
     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);
+    KelvinMatrix const D = degradation * C_tensile + C_compressive;
+    return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
+                           elastic_energy);
 }
 
 }  // namespace Phasefield

ProcessLib/PhaseField/PhaseFieldFEM-impl.h

@@ -97,8 +97,7 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
             _process_data.energy_split_model);
 
         auto& sigma = _ip_data[ip].sigma;
-        auto const& C_tensile = _ip_data[ip].C_tensile;
-        auto const& C_compressive = _ip_data[ip].C_compressive;
+        auto const& D = _ip_data[ip].D;
 
         typename ShapeMatricesType::template MatrixType<DisplacementDim,
                                                         displacement_size>
@@ -120,8 +119,7 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
             (B.transpose() * sigma - N_u.transpose() * rho_sr * b -
              local_pressure * N_u.transpose() * dNdx * d) *
             w;
-        local_Jac.noalias() +=
-            B.transpose() * (degradation * C_tensile + C_compressive) * B * w;
+        local_Jac.noalias() += B.transpose() * D * B * w;
     }
 }
 

ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h

@@ -10,10 +10,9 @@
  */
 #pragma once
 
-#include "ThermoMechanicalPhaseFieldFEM.h"
-
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
+#include "ThermoMechanicalPhaseFieldFEM.h"
 
 namespace ProcessLib
 {
@@ -93,9 +92,7 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
 
         auto& eps = _ip_data[ip].eps;
 
-        auto const& C_tensile = _ip_data[ip].C_tensile;
-        auto const& C_compressive = _ip_data[ip].C_compressive;
-
+        auto const& D = _ip_data[ip].D;
         auto const& sigma = _ip_data[ip].sigma;
 
         auto const k = _process_data.residual_stiffness(t, x_position)[0];
@@ -113,7 +110,6 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
         double const d_ip = N.dot(d);
         double const degradation = d_ip * d_ip * (1 - k) + k;
 
-        auto const C_eff = degradation * C_tensile + C_compressive;
         eps.noalias() = B * u;
         _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u, alpha,
                                                 delta_T, degradation);
@@ -131,7 +127,7 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
                 .noalias() = N;
         }
 
-        local_Jac.noalias() += B.transpose() * C_eff * B * w;
+        local_Jac.noalias() += B.transpose() * D * B * w;
 
         local_rhs.noalias() -=
             (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;

ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h

@@ -60,7 +60,7 @@ struct IntegrationPointData final
         DisplacementDim>::MaterialStateVariables>
         material_state_variables;
 
-    typename BMatricesType::KelvinMatrixType C_tensile, C_compressive;
+    typename BMatricesType::KelvinMatrixType D;
     double integration_weight;
 
     void pushBackState()
@@ -91,8 +91,8 @@ 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) =
+        std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
+                 elastic_energy) =
             MaterialLib::Solids::Phasefield::calculateVolDevDegradedStress<
                 DisplacementDim>(degradation, bulk_modulus, mu, eps_m);
     }
@@ -206,9 +206,8 @@ public:
             ip_data.eps.setZero(kelvin_vector_size);
             ip_data.eps_prev.resize(kelvin_vector_size);
             ip_data.eps_m.setZero(kelvin_vector_size);
-            ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
-            ip_data.C_compressive.setZero(kelvin_vector_size,
-                                          kelvin_vector_size);
+
+            ip_data.D.setZero(kelvin_vector_size, kelvin_vector_size);
             ip_data.sigma_tensile.setZero(kelvin_vector_size);
             ip_data.sigma_compressive.setZero(kelvin_vector_size);
             ip_data.heatflux.setZero(DisplacementDim);