Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <limits>
       8             : #include <memory>
       9             : #include <pup.h>
      10             : 
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : #include "Options/String.hpp"
      13             : #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp"
      14             : #include "Utilities/Gsl.hpp"
      15             : #include "Utilities/Serialization/CharmPupable.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : 
      18             : /// \cond
      19             : class DataVector;
      20             : /// \endcond
      21             : 
      22             : namespace Elasticity {
      23             : namespace ConstitutiveRelations {
      24             : 
      25             : /*!
      26             :  * \brief An isotropic and homogeneous material
      27             :  *
      28             :  * \details For an isotropic and homogeneous material the linear constitutive
      29             :  * relation \f$T^{ij}=-Y^{ijkl}S_{kl}\f$ reduces to
      30             :  *
      31             :  * \f[
      32             :  * Y^{ijkl} = \lambda \delta^{ij}\delta^{kl} + \mu
      33             :  * \left(\delta^{ik}\delta^{jl} + \delta^{il}\delta^{jk}\right) \\
      34             :  * \implies \quad T^{ij} = -\lambda \mathrm{Tr}(S) \delta^{ij} - 2\mu S^{ij}
      35             :  * \f]
      36             :  *
      37             :  * with the _Lamé parameter_ \f$\lambda\f$ and the _shear modulus_ (or
      38             :  * _rigidity_) \f$\mu\f$. In the parametrization chosen in this implementation
      39             :  * we use the _bulk modulus_ (or _incompressibility_)
      40             :  *
      41             :  * \f[
      42             :  * K=\lambda + \frac{2}{3}\mu
      43             :  * \f]
      44             :  *
      45             :  * instead of the Lamé parameter. In this parametrization the
      46             :  * stress-strain relation
      47             :  *
      48             :  * \f[
      49             :  * T^{ij} = -K \mathrm{Tr}(S) \delta^{ij} - 2\mu\left(S^{ij} -
      50             :  * \frac{1}{3}\mathrm{Tr}(S)\delta^{ij}\right)
      51             :  * \f]
      52             :  *
      53             :  * decomposes into a scalar and a traceless part (Eq. 11.18 in
      54             :  * \cite ThorneBlandford2017). Parameters also often used in this context are
      55             :  * the _Young's modulus_
      56             :  *
      57             :  * \f[
      58             :  * E=\frac{9K\mu}{3K+\mu}=\frac{\mu(3\lambda+2\mu)}{\lambda+\mu}
      59             :  * \f]
      60             :  *
      61             :  * and the _Poisson ratio_
      62             :  *
      63             :  * \f[
      64             :  * \nu=\frac{3K-2\mu}{2(3K+\mu)}=\frac{\lambda}{2(\lambda+\mu)}\text{.}
      65             :  * \f]
      66             :  *
      67             :  * Inversely, these relations read:
      68             :  *
      69             :  * \f[
      70             :  * K =\frac{E}{3(1-2\nu)}, \quad
      71             :  * \lambda =\frac{E\nu}{(1+\nu)(1-2\nu)}, \quad
      72             :  * \mu =\frac{E}{2(1+\nu)}
      73             :  * \f]
      74             :  *
      75             :  * **In two dimensions** this implementation reduces to the plane-stress
      76             :  * approximation. We assume that all stresses apply in the plane of the
      77             :  * computational domain, which corresponds to scenarios of in-plane stretching
      78             :  * and shearing of thin slabs of material. Since orthogonal stresses vanish as
      79             :  * \f$T^{i3}=0=T^{3i}\f$ we find \f$\mathrm{Tr}(S)=\frac{2\mu}{\lambda +
      80             :  * 2\mu}\mathrm{Tr}^{(2)}(S)\f$, where \f$\mathrm{Tr}^{(2)}\f$ denotes that the
      81             :  * trace only applies to the two dimensions within the plane. The constitutive
      82             :  * relation thus reduces to
      83             :  *
      84             :  * \f{align}
      85             :  * T^{ij}=&-\frac{2\lambda\mu}{\lambda + 2\mu}\mathrm{Tr}^{(2)}(S)\delta^{ij} -
      86             :  * 2\mu S^{ij} \\
      87             :  * =&-\frac{E\nu}{1-\nu^2}\mathrm{Tr}^{(2)}(S)\delta^{ij} -
      88             :  * \frac{E}{1+\nu}S^{ij}
      89             :  * \f}
      90             :  *
      91             :  * which is non-zero only in the directions of the plane. Since the stresses
      92             :  * are also assumed to be constant along the thickness of the plane
      93             :  * \f$\partial_3T^{ij}=0\f$ the elasticity problem \f$-\partial_i T^{ij}=F^j\f$
      94             :  * reduces to two dimensions.
      95             :  */
      96             : template <size_t Dim>
      97           1 : class IsotropicHomogeneous : public ConstitutiveRelation<Dim> {
      98             :  public:
      99           0 :   static constexpr size_t volume_dim = Dim;
     100             : 
     101           0 :   struct BulkModulus {
     102           0 :     using type = double;
     103           0 :     static constexpr Options::String help = {
     104             :         "The incompressibility of the material"};
     105           0 :     static type lower_bound() { return 0.0; }
     106             :   };
     107             : 
     108           0 :   struct ShearModulus {
     109           0 :     using type = double;
     110           0 :     static constexpr Options::String help = {"The rigidity of the material"};
     111           0 :     static type lower_bound() { return 0.0; }
     112             :   };
     113             : 
     114           0 :   using options = tmpl::list<BulkModulus, ShearModulus>;
     115             : 
     116           0 :   static constexpr Options::String help = {
     117             :       "A constitutive relation that describes an isotropic, homogeneous "
     118             :       "material in terms of two elastic moduli. These bulk and shear moduli "
     119             :       "indicate the material's resistance to volume and shape changes, "
     120             :       "respectively. Both are measured in units of stress, typically Pascals."};
     121             : 
     122           0 :   IsotropicHomogeneous() = default;
     123           0 :   IsotropicHomogeneous(const IsotropicHomogeneous&) = default;
     124           0 :   IsotropicHomogeneous& operator=(const IsotropicHomogeneous&) = default;
     125           0 :   IsotropicHomogeneous(IsotropicHomogeneous&&) = default;
     126           0 :   IsotropicHomogeneous& operator=(IsotropicHomogeneous&&) = default;
     127           0 :   ~IsotropicHomogeneous() override = default;
     128             : 
     129           0 :   IsotropicHomogeneous(double bulk_modulus, double shear_modulus);
     130             : 
     131           1 :   std::unique_ptr<ConstitutiveRelation<Dim>> get_clone() const override;
     132             : 
     133             :   /// The constitutive relation that characterizes the elastic properties of a
     134             :   /// material
     135           1 :   void stress(gsl::not_null<tnsr::II<DataVector, Dim>*> stress,
     136             :               const tnsr::ii<DataVector, Dim>& strain,
     137             :               const tnsr::I<DataVector, Dim>& x) const override;
     138             : 
     139             :   /// The bulk modulus (or incompressibility) \f$K\f$
     140           1 :   double bulk_modulus() const;
     141             :   /// The shear modulus (or rigidity) \f$\mu\f$
     142           1 :   double shear_modulus() const;
     143             :   /// The Lamé parameter \f$\lambda\f$
     144           1 :   double lame_parameter() const;
     145             :   /// The Young's modulus \f$E\f$
     146           1 :   double youngs_modulus() const;
     147             :   /// The Poisson ratio \f$\nu\f$
     148           1 :   double poisson_ratio() const;
     149             : 
     150             :   // clang-tidy: no runtime references
     151           0 :   void pup(PUP::er& /*p*/) override;  //  NOLINT
     152             : 
     153           0 :   explicit IsotropicHomogeneous(CkMigrateMessage* /*unused*/) {}
     154             : 
     155           0 :   WRAPPED_PUPable_decl_base_template(  // NOLINT
     156             :       SINGLE_ARG(ConstitutiveRelation<Dim>), IsotropicHomogeneous);
     157             : 
     158             :  private:
     159           0 :   double bulk_modulus_ = std::numeric_limits<double>::signaling_NaN();
     160           0 :   double shear_modulus_ = std::numeric_limits<double>::signaling_NaN();
     161             : };
     162             : 
     163             : template <size_t Dim>
     164           0 : bool operator==(const IsotropicHomogeneous<Dim>& lhs,
     165             :                 const IsotropicHomogeneous<Dim>& rhs);
     166             : template <size_t Dim>
     167           0 : bool operator!=(const IsotropicHomogeneous<Dim>& lhs,
     168             :                 const IsotropicHomogeneous<Dim>& rhs);
     169             : 
     170             : /// \cond
     171             : template <size_t Dim>
     172             : PUP::able::PUP_ID IsotropicHomogeneous<Dim>::my_PUP_ID = 0;
     173             : /// \endcond
     174             : 
     175             : }  // namespace ConstitutiveRelations
     176             : }  // namespace Elasticity