arrudaboyce.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <ikarus/finiteelements/mechanics/materials/hyperelastic/deviatoric/deviatoricinvariants.hh>
#include <ikarus/finiteelements/mechanics/materials/materialhelpers.hh>
#include <ikarus/finiteelements/mechanics/materials/tags.hh>
#include <ikarus/utils/tensorutils.hh>
 
namespace Ikarus {
 
struct ArrudaBoyceMatParameters
{
  double mu;      
  double lambdaM; 
};
} // namespace Ikarus
 
namespace Ikarus::Materials {
 
template <typename ST_>
struct ArrudaBoyceT
{
  static constexpr int dim = 3;
  using ScalarType         = ST_;
 
  template <typename ST = ScalarType>
  using PrincipalStretches = Eigen::Vector<ST, 3>;
  template <typename ST = ScalarType>
  using Invariants = PrincipalStretches<ST>;
 
  template <typename ST = ScalarType>
  using FirstDerivative = Eigen::Vector<ST, dim>;
  template <typename ST = ScalarType>
  using SecondDerivative = Eigen::Matrix<ST, dim, dim>;
 
  using MaterialParameters      = ArrudaBoyceMatParameters;
  static constexpr int numTerms = 5;
 
  static constexpr auto stretchTag = PrincipalStretchTags::deviatoric;
 
  [[nodiscard]] constexpr static std::string name() noexcept { return "ArrudaBoyce"; }
 
  explicit ArrudaBoyceT(const MaterialParameters& matPar)
      : matPar_{matPar} {}
 
  MaterialParameters materialParametersImpl() const { return matPar_; }
 
  template <typename ST>
  ST storedEnergyImpl(const PrincipalStretches<ST>& lambda) const {
    const Invariants<ST>& invariants = Impl::invariants(lambda);
    ST energy{0.0};
    const auto& devInvariants = DeviatoricInvariants<PrincipalStretches<ST>>(lambda);
    const auto W1             = devInvariants.value()[0];
    const auto mu             = matPar_.mu;
    const auto lambdaM        = matPar_.lambdaM;
    const auto beta           = 1 / pow(lambdaM, 2.0);
 
    for (const auto i : parameterRange())
      energy += alphas_[i] * pow(beta, i) * (pow(W1, i + 1) - pow(3, i + 1));
    energy *= mu;
 
    return energy;
  }
 
  template <typename ST>
  FirstDerivative<ST> firstDerivativeImpl(const PrincipalStretches<ST>& lambda) const {
    const Invariants<ST>& invariants = Impl::invariants(lambda);
    auto dWdLambda                   = FirstDerivative<ST>::Zero().eval();
 
    const auto& devInvariants = DeviatoricInvariants<PrincipalStretches<ST>>(lambda);
    const auto W1             = devInvariants.value()[0];
    const auto& dW1dLambda    = devInvariants.firstDerivative().first;
    const auto mu             = matPar_.mu;
    const auto lambdaM        = matPar_.lambdaM;
    const auto beta           = 1 / pow(lambdaM, 2.0);
 
    for (const auto j : parameterRange())
      for (const auto k : dimensionRange())
        dWdLambda[k] += mu * alphas_[j] * pow(beta, j) * pow(W1, j) * dW1dLambda[k] * (j + 1);
 
    return dWdLambda;
  }
 
  template <typename ST>
  SecondDerivative<ST> secondDerivativeImpl(const PrincipalStretches<ST>& lambda) const {
    const Invariants<ST>& invariants = Impl::invariants(lambda);
    auto dS                          = SecondDerivative<ST>::Zero().eval();
 
    const auto& devInvariants = DeviatoricInvariants<PrincipalStretches<ST>>(lambda);
    const auto W1             = devInvariants.value()[0];
    const auto& dW1dLambda    = devInvariants.firstDerivative().first;
    const auto& ddW1dLambda   = devInvariants.secondDerivative().first;
    const auto mu             = matPar_.mu;
    const auto lambdaM        = matPar_.lambdaM;
    const auto beta           = 1 / pow(lambdaM, 2.0);
 
    const auto dW1dLambdaDyad = dyadic(dW1dLambda, dW1dLambda);
    for (const auto p : parameterRange()) {
      const auto factor1 = mu * alphas_[p] * pow(beta, p);
      const auto factor2 = pow(W1, p) * ddW1dLambda * (p + 1);
      const auto factor3 = pow(W1, p - 1) * dW1dLambdaDyad * p * (p + 1);
      dS += factor1 * (factor2 + factor3);
      dS.diagonal() -= (dW1dLambda.array() / lambda.array() * factor1 * pow(W1, p) * (p + 1)).matrix();
    }
 
    return dS;
  }
 
  template <typename STO>
  auto rebind() const {
    return ArrudaBoyceT<STO>(matPar_);
  }
 
private:
  MaterialParameters matPar_;
  static constexpr std::array<double, numTerms> alphas_ = {0.5, 1.0 / 20.0, 11.0 / 1050.0, 19.0 / 7000.0,
                                                           519.0 / 673750.0};
 
  inline static constexpr auto parameterRange() { return Dune::range(numTerms); }
  inline static constexpr auto dimensionRange() { return Dune::range(dim); }
};
 
using ArrudaBoyce = ArrudaBoyceT<double>;
 
} // namespace Ikarus::Materials