[PL] SD; Non-equilibrium initial stress state.

ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp

@@ -101,9 +101,25 @@ std::unique_ptr<Process> createSmallDeformationProcess(
         config.getConfigParameter<double>(
             "reference_temperature", std::numeric_limits<double>::quiet_NaN());
 
+    // Non-equilibrium variables
+    ParameterLib::Parameter<double> const* nonequilibrium_stress = nullptr;
+    const auto& nonequilibrium_state_variables_config =
+        //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION__nonequilibrium_state_variables}
+        config.getConfigSubtreeOptional("nonequilibrium_state_variables");
+    if (nonequilibrium_state_variables_config)
+    {
+        nonequilibrium_stress = &ParameterLib::findParameter<double>(
+            *nonequilibrium_state_variables_config,
+            //! \ogs_file_param_special{prj__processes__process__SMALL_DEFORMATION__nonequilibrium_state_variables__stress}
+            "stress", parameters,
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value);
+    }
+
     SmallDeformationProcessData<DisplacementDim> process_data{
-        materialIDs(mesh), std::move(solid_constitutive_relations),
-        solid_density, specific_body_force, reference_temperature};
+        materialIDs(mesh),     std::move(solid_constitutive_relations),
+        solid_density,         specific_body_force,
+        reference_temperature, nonequilibrium_stress};
 
     SecondaryVariableCollection secondary_variables;
 

ProcessLib/SmallDeformation/SmallDeformationFEM.h

@@ -48,6 +48,13 @@ struct IntegrationPointData final
 
     typename BMatricesType::KelvinVectorType sigma, sigma_prev;
     typename BMatricesType::KelvinVectorType eps, eps_prev;
+
+    /// Non-equilibrium initial stress.
+    typename BMatricesType::KelvinVectorType sigma_neq =
+        BMatricesType::KelvinVectorType::Zero(
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value);
+
     double free_energy_density = 0;
 
     MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
@@ -132,6 +139,8 @@ public:
                 _process_data.material_ids,
                 e.getID());
 
+        ParameterLib::SpatialPosition x_position;
+
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             _ip_data.emplace_back(solid_material);
@@ -147,8 +156,30 @@ public:
             static const int kelvin_vector_size =
                 MathLib::KelvinVector::KelvinVectorDimensions<
                     DisplacementDim>::value;
+
             // Initialize current time step values
-            ip_data.sigma.setZero(kelvin_vector_size);
+            if (_process_data.nonequilibrium_stress)
+            {
+                // Computation of non-equilibrium stress.
+                x_position.setCoordinates(MathLib::Point3d(
+                    interpolateCoordinates<ShapeFunction, ShapeMatricesType>(
+                        e, ip_data.N)));
+                std::vector<double> sigma_neq_data =
+                    (*_process_data.nonequilibrium_stress)(
+                        std::numeric_limits<
+                            double>::quiet_NaN() /* time independent */,
+                        x_position);
+                ip_data.sigma_neq =
+                    Eigen::Map<typename BMatricesType::KelvinVectorType>(
+                        sigma_neq_data.data(),
+                        MathLib::KelvinVector::KelvinVectorDimensions<
+                            DisplacementDim>::value,
+                        1);
+            }
+            // Initialization from non-equilibrium sigma, which is zero by
+            // default, or is set to some value.
+            ip_data.sigma = ip_data.sigma_neq;
+
             ip_data.eps.setZero(kelvin_vector_size);
 
             // Previous time step values are not initialized and are set later.
@@ -246,6 +277,7 @@ public:
 
             auto const& eps_prev = _ip_data[ip].eps_prev;
             auto const& sigma_prev = _ip_data[ip].sigma_prev;
+            auto const& sigma_neq = _ip_data[ip].sigma_neq;
 
             auto& eps = _ip_data[ip].eps;
             auto& sigma = _ip_data[ip].sigma;
@@ -270,8 +302,9 @@ public:
 
             auto const rho = _process_data.solid_density(t, x_position)[0];
             auto const& b = _process_data.specific_body_force;
-            local_b.noalias() -=
-                (B.transpose() * sigma - N_u_op.transpose() * rho * b) * w;
+            local_b.noalias() -= (B.transpose() * (sigma - sigma_neq) -
+                                  N_u_op.transpose() * rho * b) *
+                                 w;
             local_Jac.noalias() += B.transpose() * C * B * w;
         }
     }

ProcessLib/SmallDeformation/SmallDeformationProcessData.h

@@ -40,10 +40,12 @@ struct SmallDeformationProcessData
         ParameterLib::Parameter<double> const& solid_density_,
         Eigen::Matrix<double, DisplacementDim, 1>
             specific_body_force_,
-        double const reference_temperature_)
+        double const reference_temperature_,
+        ParameterLib::Parameter<double> const* const nonequilibrium_stress_)
         : material_ids(material_ids_),
           solid_materials{std::move(solid_materials_)},
           solid_density(solid_density_),
+          nonequilibrium_stress(nonequilibrium_stress_),
           specific_body_force(std::move(specific_body_force_)),
           reference_temperature(reference_temperature_)
     {
@@ -68,6 +70,8 @@ struct SmallDeformationProcessData
         solid_materials;
     /// Solid's density. A scalar quantity, ParameterLib::Parameter<double>.
     ParameterLib::Parameter<double> const& solid_density;
+
+    ParameterLib::Parameter<double> const* const nonequilibrium_stress;
     /// Specific body forces applied to the solid.
     /// It is usually used to apply gravitational forces.
     /// A vector of displacement dimension's length.