solid_material.hpp

Go to the documentation of this file.

// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
// other Serac Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include "serac/numerics/functional/functional.hpp"
 
namespace serac::solid_mechanics {
 
struct LinearIsotropic {
  using State = Empty;  
 
  template <typename T, int dim>
  SERAC_HOST_DEVICE auto operator()(State& /* state */, const tensor<T, dim, dim>& du_dX) const
  {
    auto I       = Identity<dim>();
    auto lambda  = K - (2.0 / 3.0) * G;
    auto epsilon = 0.5 * (transpose(du_dX) + du_dX);
    return lambda * tr(epsilon) * I + 2.0 * G * epsilon;
  }
 
  double density;  
  double K;        
  double G;        
};
 
template <typename T, int dim>
auto greenStrain(const tensor<T, dim, dim>& grad_u)
{
  return 0.5 * (grad_u + transpose(grad_u) + dot(transpose(grad_u), grad_u));
}
 
struct StVenantKirchhoff {
  using State = Empty;  
 
  template <typename T, int dim>
  auto operator()(State&, const tensor<T, dim, dim>& grad_u) const
  {
    static constexpr auto I = Identity<dim>();
    auto                  F = grad_u + I;
    const auto            E = greenStrain(grad_u);
 
    // stress
    const auto S     = K * tr(E) * I + 2.0 * G * dev(E);
    const auto P     = dot(F, S);
    const auto sigma = dot(P, transpose(F)) / det(F);
 
    return sigma;
  }
 
  double density;  
  double K;        
  double G;        
};
 
struct NeoHookean {
  using State = Empty;  
 
  template <typename T, int dim>
  SERAC_HOST_DEVICE auto operator()(State& /* state */, const tensor<T, dim, dim>& du_dX) const
  {
    using std::log1p;
    constexpr auto I         = Identity<dim>();
    auto           lambda    = K - (2.0 / 3.0) * G;
    auto           B_minus_I = du_dX * transpose(du_dX) + transpose(du_dX) + du_dX;
    auto           J_minus_1 = detApIm1(du_dX);
    auto           J         = J_minus_1 + 1;
    return (lambda * log1p(J_minus_1) * I + G * B_minus_I) / J;
  }
 
  double density;  
  double K;        
  double G;        
};
 
struct PowerLawHardening {
  double sigma_y;  
  double n;        
  double eps0;     
 
  template <typename T>
  auto operator()(const T accumulated_plastic_strain) const
  {
    using std::pow;
    return sigma_y * pow(1.0 + accumulated_plastic_strain / eps0, 1.0 / n);
  };
};
 
struct VoceHardening {
  double sigma_y;          
  double sigma_sat;        
  double strain_constant;  
 
  template <typename T>
  auto operator()(const T accumulated_plastic_strain) const
  {
    using std::exp;
    return sigma_sat - (sigma_sat - sigma_y) * exp(-accumulated_plastic_strain / strain_constant);
  };
};
 
template <typename HardeningType>
struct J2Nonlinear {
  static constexpr int    dim = 3;      
  static constexpr double tol = 1e-10;  
 
  double        E;          
  double        nu;         
  HardeningType hardening;  
  double        density;    
 
  struct State {
    tensor<double, dim, dim> plastic_strain;              
    double                   accumulated_plastic_strain;  
  };
 
  template <typename T>
  auto operator()(State& state, const T du_dX) const
  {
    using std::sqrt;
    constexpr auto I = Identity<dim>();
    const double   K = E / (3.0 * (1.0 - 2.0 * nu));
    const double   G = 0.5 * E / (1.0 + nu);
 
    // (i) elastic predictor
    auto el_strain = sym(du_dX) - state.plastic_strain;
    auto p         = K * tr(el_strain);
    auto s         = 2.0 * G * dev(el_strain);
    auto q         = sqrt(1.5) * norm(s);
 
    // (ii) admissibility
    const double eqps_old = state.accumulated_plastic_strain;
    auto         residual = [eqps_old, G, *this](auto delta_eqps, auto trial_mises) {
      return trial_mises - 3.0 * G * delta_eqps - this->hardening(eqps_old + delta_eqps);
    };
    if (residual(0.0, get_value(q)) > tol * hardening.sigma_y) {
      // (iii) return mapping
 
      // Note the tolerance for convergence is the same as the tolerance for entering the return map.
      // This ensures that if the constitutive update is called again with the updated internal
      // variables, the return map won't be repeated.
      ScalarSolverOptions opts{.xtol = 0, .rtol = tol * hardening.sigma_y, .max_iter = 25};
      double              lower_bound = 0.0;
      double              upper_bound = (get_value(q) - hardening(eqps_old)) / (3.0 * G);
      auto [delta_eqps, status]       = solve_scalar_equation(residual, 0.0, lower_bound, upper_bound, opts, q);
 
      auto Np = 1.5 * s / q;
 
      s = s - 2.0 * G * delta_eqps * Np;
      state.accumulated_plastic_strain += get_value(delta_eqps);
      state.plastic_strain += get_value(delta_eqps) * get_value(Np);
    }
 
    return s + p * I;
  }
};
 
struct J2 {
  static constexpr int dim = 3;
 
  double E;        
  double nu;       
  double Hi;       
  double Hk;       
  double sigma_y;  
  double density;  
 
  struct State {
    tensor<double, dim, dim> beta;                        
    tensor<double, dim, dim> plastic_strain;              
    double                   accumulated_plastic_strain;  
  };
 
  template <typename T>
  auto operator()(State& state, const T du_dX) const
  {
    using std::sqrt;
    constexpr auto I = Identity<3>();
    const double   K = E / (3.0 * (1.0 - 2.0 * nu));
    const double   G = 0.5 * E / (1.0 + nu);
 
    //
    // see pg. 260, box 7.5,
    // in "Computational Methods for Plasticity"
    //
 
    // (i) elastic predictor
    auto el_strain = sym(du_dX) - state.plastic_strain;
    auto p         = K * tr(el_strain);
    auto s         = 2.0 * G * dev(el_strain);
    auto eta       = s - state.beta;
    auto q         = sqrt(3.0 / 2.0) * norm(eta);
    auto phi       = q - (sigma_y + Hi * state.accumulated_plastic_strain);
 
    // (ii) admissibility
    if (phi > 0.0) {
      // see (7.207) on pg. 261
      auto plastic_strain_inc = phi / (3 * G + Hk + Hi);
 
      // from here on, only normalize(eta) is required
      // so we overwrite eta with its normalized version
      eta = normalize(eta);
 
      // (iii) return mapping
      s = s - sqrt(6.0) * G * plastic_strain_inc * eta;
      state.accumulated_plastic_strain += get_value(plastic_strain_inc);
      state.plastic_strain += sqrt(3.0 / 2.0) * get_value(plastic_strain_inc) * get_value(eta);
      state.beta = state.beta + sqrt(2.0 / 3.0) * Hk * get_value(plastic_strain_inc) * get_value(eta);
    }
 
    return s + p * I;
  }
};
 
template <typename T1, typename T2, int dim>
auto KirchhoffToPiola(const tensor<T1, dim, dim>& kirchhoff_stress, const tensor<T2, dim, dim>& displacement_gradient)
{
  return transpose(dot(inv(displacement_gradient + Identity<dim>()), kirchhoff_stress));
}
 
template <typename T1, typename T2, int dim>
auto CauchyToPiola(const tensor<T1, dim, dim>& cauchy_stress, const tensor<T2, dim, dim>& displacement_gradient)
{
  auto kirchhoff_stress = det(displacement_gradient + Identity<dim>()) * cauchy_stress;
  return KirchhoffToPiola(kirchhoff_stress, displacement_gradient);
}
 
template <int dim>
struct ConstantBodyForce {
  tensor<double, dim> force_;
 
  template <typename T>
  SERAC_HOST_DEVICE tensor<double, dim> operator()(const tensor<T, dim>& /* x */, const double /* t */) const
  {
    return force_;
  }
};
 
template <int dim>
struct ConstantTraction {
  tensor<double, dim> traction_;
 
  template <typename T>
  SERAC_HOST_DEVICE tensor<double, dim> operator()(const tensor<T, dim>& /* x */, const tensor<T, dim>& /* n */,
                                                   const double /* t */) const
  {
    return traction_;
  }
};
 
}  // namespace serac::solid_mechanics