LinearElasticTransverseIsotropic.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include "MechanicsBase.h"
#include "ParameterLib/Parameter.h"
 
namespace MaterialLib
{
namespace Solids
{
template <int DisplacementDim>
class LinearElasticTransverseIsotropic : public MechanicsBase<DisplacementDim>
{
public:
    static int const KelvinVectorSize =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
    using KelvinVector =
        MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
    using KelvinMatrix =
        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
 
    using P = ParameterLib::Parameter<double>;
 
    LinearElasticTransverseIsotropic(
        P const& E_i, P const& E_a, P const& nu_ii, P const& nu_ia,
        P const& G_ia,
        std::optional<ParameterLib::CoordinateSystem> const&
            local_coordinate_system)
        : E_i_p_(E_i),
          E_a_p_(E_a),
          nu_ii_p_(nu_ii),
          nu_ia_p_(nu_ia),
          G_ia_p_(G_ia),
          local_coordinate_system_(local_coordinate_system)
    {
    }
 
    double computeFreeEnergyDensity(
        double const /*t*/,
        ParameterLib::SpatialPosition const& /*x*/,
        double const /*dt*/,
        KelvinVector const& eps,
        KelvinVector const& sigma,
        typename MechanicsBase<DisplacementDim>::
            MaterialStateVariables const& /* material_state_variables */)
        const override
    {
        return eps.dot(sigma) / 2;
    }
 
    std::optional<
        std::tuple<typename MechanicsBase<DisplacementDim>::KelvinVector,
                   std::unique_ptr<typename MechanicsBase<
                       DisplacementDim>::MaterialStateVariables>,
                   typename MechanicsBase<DisplacementDim>::KelvinMatrix>>
    integrateStress(
        MaterialPropertyLib::VariableArray const& variable_array_prev,
        MaterialPropertyLib::VariableArray const& variable_array,
        double const t, ParameterLib::SpatialPosition const& x, double const dt,
        typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
            material_state_variables) const override;
 
    KelvinMatrix getElasticTensor(double const t,
                                  ParameterLib::SpatialPosition const& x,
                                  double const T) const;
 
    double getBulkModulus(double const t,
                          ParameterLib::SpatialPosition const& x,
                          KelvinMatrix const* const /*C*/) const override
    {
        const double E_i = E_i_p_(t, x)[0];
        const double E_a = E_a_p_(t, x)[0];
        const double nu_i = nu_ii_p_(t, x)[0];
        const double nu_ia = nu_ia_p_(t, x)[0];
 
        // Average Young's modulus
        double const E_av = 2. * E_i / 3. + E_a / 3.;
 
        // Poisson ratio in the plane of isotropy, due to the strain
        // perpendicular to the plane of isotropy.  nu_ai=nu_ia*E_a/E_i
        const double nu_ai = nu_ia * E_a / E_i;
 
        // Average Poisson ratio
        //    12     13    21   23   31    32
        //    ai     ai    ia   ii   ia    ii
        const double nu_av = (nu_ai + nu_ia + nu_i) / 3.0;
 
        return E_av / 3 / (1 - 2 * nu_av);
    }
 
protected:
    P const& E_i_p_;
    P const& E_a_p_;
    P const& nu_ii_p_;
    P const& nu_ia_p_;
    P const& G_ia_p_;
    std::optional<ParameterLib::CoordinateSystem> const&
        local_coordinate_system_;
 
private:
    KelvinMatrix getElasticTensorLeftTopCorner(
        double const t, ParameterLib::SpatialPosition const& x) const;
};
 
}  // namespace Solids
}  // namespace MaterialLib