MenetreyWillam.h

Go to the documentation of this file.

/* ---------------------------------------------------------------------
 *                                       _
 *  _ __ ___   __ _ _ __ _ __ ___   ___ | |_
 * | '_ ` _ \ / _` | '__| '_ ` _ \ / _ \| __|
 * | | | | | | (_| | |  | | | | | | (_) | |_
 * |_| |_| |_|\__,_|_|  |_| |_| |_|\___/ \__|
 *
 * Unit of Strength of Materials and Structural Analysis
 * University of Innsbruck,
 * 2020 - today
 *
 * festigkeitslehre@uibk.ac.at
 *
 * Matthias Neuner matthias.neuner@uibk.ac.at
 * Magdalena Schreter magdalena.schreter@uibk.ac.at
 *
 * This file is part of the MAteRialMOdellingToolbox (marmot).
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of marmot.
 * ---------------------------------------------------------------------
 */
 
#pragma once
#include "Marmot/HaighWestergaard.h"
#include <utility>
 
namespace Marmot {
  namespace ContinuumMechanics::CommonConstitutiveModels {
 
    // clang-format off
    // clang-format on
    class MenetreyWillam {
    public:
      struct MenetreyWillamParameters {
        double Af; 
        double Bf; 
        double Cf; 
        double m;  
        double e;  
      } param;
 
      // Reduction of the generalized failure criterion to a specific type.
      enum class MenetreyWillamType {
        Mises,         
        Rankine,       
        DruckerPrager, 
        MohrCoulomb    
      };
 
      MenetreyWillam( const double              ft,
                      const MenetreyWillamType& type = MenetreyWillamType::Rankine,
                      const double              fc   = 0 );
 
      void setParameters( const double ft, const double fc, const MenetreyWillamType& type );
 
      template < typename T >
      T polarRadius( const double& theta ) const
      {
        return polarRadius( theta, param.e );
      }
 
      template < typename T >
      static T polarRadius( const T& theta, const double& e )
      {
        // computes the deviatoric shape roundness for a given eccentricity (e) at a certain
        // position (theta) the numerator and denominator are stored for performance reasons, as
        // they are also needed for the derivative dRdTheta
        if ( e >= 1.0 )
          return 1;
 
        const double e2        = e * e;
        const T      cosTheta  = cos( theta );
        const T      cos2Theta = cosTheta * cosTheta;
 
        const T numerator   = ( 4. * ( 1. - e2 ) * cos2Theta + ( 2. * e - 1. ) * ( 2. * e - 1. ) );
        const T denominator = ( 2. * ( 1. - e2 ) * cosTheta +
                                ( 2. * e - 1. ) * sqrt( 4. * ( 1. - e2 ) * cos2Theta + 5. * e2 - 4. * e ) );
 
        const T r = numerator / denominator;
 
        return r;
      }
      template < typename T >
      std::pair< T, T > dPolarRadius_dTheta( const T& theta ) const
      {
        return dPolarRadius_dTheta( theta, param.e );
      }
 
      template < typename T >
      static std::pair< T, T > dPolarRadius_dTheta( const T& theta, const double& e )
      {
        //
        // computes the deviatoric shape roundness for a given eccentricity (e) at a certain
        // position (theta) the numerator and denominator are stored for performance reasons, as
        // they are also needed for the derivative dRdTheta
        T r, dRdTheta;
 
        if ( e >= 1.0 ) {
          r        = 1;
          dRdTheta = 0;
          return std::make_pair( r, dRdTheta );
        }
 
        const double e2        = e * e;
        const T      cosTheta  = cos( theta );
        const T      cos2Theta = cosTheta * cosTheta;
 
        T numerator   = ( 4. * ( 1. - e2 ) * cos2Theta + ( 2. * e - 1. ) * ( 2. * e - 1. ) );
        T denominator = ( 2. * ( 1. - e2 ) * cosTheta +
                          ( 2. * e - 1. ) * sqrt( 4. * ( 1. - e2 ) * cos2Theta + 5. * e2 - 4. * e ) );
 
        r = numerator / denominator;
 
        // compute the derivative
        const T sinTheta = sin( theta );
 
        const T a          = numerator;
        const T b          = 1. / denominator;
        const T aux        = 4. * ( 1. - e2 ) * cos2Theta + 5. * e2 - 4. * e;
        const T dAuxdTheta = 4. * ( 1. - e2 ) * 2. * cosTheta * -sinTheta;
        const T dadTheta   = 2. * cosTheta * ( -sinTheta ) * 4. * ( 1. - e2 );
        const T dbdTheta   = -pow( denominator, -2. ) * ( 2. * ( 1. - e2 ) * -sinTheta +
                                                        ( 2. * e - 1 ) * 1. / 2 * pow( aux, -1. / 2 ) * dAuxdTheta );
        dRdTheta           = b * dadTheta + a * dbdTheta;
 
        return { r, dRdTheta };
      }
 
      template < typename T >
      static std::tuple< T, T, T > d2PolarRadius_dTheta2( const T& theta, const double& e )
      {
 
        const auto [r, dr_dTheta] = dPolarRadius_dTheta< T >( theta, e );
        if ( e >= 1.0 ) {
          return { r, dr_dTheta, 0 };
        }
 
        const double e2   = pow( e, 2 );
        const double e2m1 = e2 - 1;
        const double tem1 = 2 * e - 1;
 
        const T cos_theta   = cos( theta );
        const T sin_theta   = sin( theta );
        const T cos_2_theta = cos_theta * cos_theta;
        const T sin_2_theta = sin_theta * sin_theta;
 
        const T x      = 5. * e2 - 4. * e - 4. * e2m1 * cos_2_theta;
        const T sqrt_x = sqrt( x );
 
        const T one_over_sqrt_x = 1. / sqrt_x;
 
        const T d2r_dTheta2 = e2m1 *
                              ( -16. * e2m1 * ( 4.0 * tem1 * one_over_sqrt_x * cos_theta + 2. ) * sin_2_theta *
                                  cos_theta / ( tem1 * sqrt_x - 2 * e2m1 * cos_theta ) -
                                8. * sin_2_theta + 8. * cos_2_theta +
                                ( pow( 2. * e - 1., 2 ) - 4 * e2m1 * cos_2_theta ) *
                                  ( 16.0 * tem1 * e2m1 * pow( x, -1.5 ) * sin_2_theta * cos_2_theta +
                                    4.0 * tem1 * one_over_sqrt_x * sin_2_theta -
                                    4.0 * tem1 * one_over_sqrt_x * cos_2_theta +
                                    e2m1 * ( 4.0 * tem1 * one_over_sqrt_x * cos_theta + 2. ) *
                                      ( 8.0 * tem1 * one_over_sqrt_x * cos_theta + 4. ) * sin_2_theta /
                                      ( tem1 * sqrt_x - 2. * e2m1 * cos_theta ) -
                                    2. * cos_theta ) /
                                  ( tem1 * sqrt_x - 2. * e2m1 * cos_theta ) ) /
                              ( tem1 * sqrt_x - 2. * e2m1 * cos_theta );
 
        return { r, dr_dTheta, d2r_dTheta2 };
      }
 
      template < typename T >
      T yieldFunction( const ContinuumMechanics::HaighWestergaard::HaighWestergaardCoordinates< T >& hw,
                       const double varEps = 0.0 ) const
      {
        const T r_ = polarRadius( hw.theta, param.e );
        if ( varEps == 0 )
          return ( param.Af * hw.rho ) * ( param.Af * hw.rho ) +
                 param.m * ( param.Bf * hw.rho * r_ + param.Cf * hw.xi ) - 1.;
        else
          return param.Af * param.Af * hw.rho * hw.rho +
                 param.m * ( std::sqrt( param.Bf * hw.rho * r_ * param.Bf * hw.rho * r_ + varEps * varEps ) +
                             param.Cf * hw.xi ) -
                 1.;
      }
 
      template < typename T >
      std::tuple< T, T, T > dYieldFunction_dHaighWestergaard(
        const ContinuumMechanics::HaighWestergaard::HaighWestergaardCoordinates< T >& hw,
        const double                                                                  varEps = 0.0 ) const
      {
        const auto [r_, dRdTheta_] = dPolarRadius_dTheta( hw.theta, param.e );
 
        T dFdXi, dFdRho, dFdTheta;
        dFdXi = param.m * param.Cf;
 
        if ( varEps == 0.0 ) {
          dFdRho   = 2. * param.Af * hw.rho * param.Af + param.m * param.Bf * r_;
          dFdTheta = param.m * param.Bf * hw.rho * dRdTheta_;
        }
        else {
          const T auxTerm1 = param.m * 0.5 *
                             pow( param.Bf * param.Bf * hw.rho * hw.rho * r_ * r_ + varEps * varEps, -0.5 ) * 2 *
                             param.Bf * hw.rho * r_ * param.Bf;
 
          dFdRho   = param.Af * param.Af * 2. * hw.rho + auxTerm1 * r_;
          dFdTheta = auxTerm1 * hw.rho * dRdTheta_;
        }
 
        return { dFdXi, dFdRho, dFdTheta };
      }
 
      static inline double abaqusMohrCoulombPotentialVarEpsToMenetreyWillam( const double varEps, const double psi )
      {
        return varEps * 2 * std::sin( psi );
      }
 
      static inline double e( const double fc, const double ft ) { return ( fc + 2 * ft ) / ( 2 * fc + ft ); }
 
      static inline double c( const double fc, const double ft )
      {
        const double phi_ = phi( fc, ft );
        return ft * ( 1 + std::sin( phi_ ) ) / ( 2 + std::cos( phi_ ) );
      }
 
      static inline double phi( const double fc, const double ft ) { return std::asin( ( fc - ft ) / ( fc + ft ) ); }
 
      static inline double ft( const double c, const double phi )
      {
        return 2 * c * std::cos( phi ) / ( 1 + std::sin( phi ) );
      }
 
      static inline double fc( const double c, const double phi )
      {
        return 2 * c * std::cos( phi ) / ( 1 - std::sin( phi ) );
      }
    };
 
  } // namespace ContinuumMechanics::CommonConstitutiveModels
} // namespace Marmot