solid_material.hpp

Go to the documentation of this file.

// Copyright (c) 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);
    return dot(F, S);
  }
 
  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 = dot(du_dX, transpose(du_dX)) + transpose(du_dX) + du_dX;
 
    auto logJ = log1p(detApIm1(du_dX));
    // Kirchoff stress, in form that avoids cancellation error when F is near I
    auto TK = lambda * logJ * I + G * B_minus_I;
 
    // Pull back to Piola
    auto F = du_dX + I;
    return dot(TK, inv(transpose(F)));
  }
 
  double density;  
  double K;        
  double G;        
};
 
struct NeoHookeanAdditiveSplit {
  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::pow;
    constexpr auto I = Identity<dim>();
    auto F = I + du_dX;
    auto J = det(F);
    auto Jm13 = pow(J, -1.0 / 3.0);
    auto F_bar = Jm13 * F;
    auto Pdev = G * Jm13 * (F_bar - 1.0 / 3.0 * inner(F_bar, F_bar) * inv(transpose(F_bar)));
 
    // any volumetric model could be substituted here
    using std::log1p;
    auto logJ = log1p(detApIm1(du_dX));
    auto Pvol = K * logJ * inv(transpose(F));
 
    return Pdev + Pvol;
  }
 
  double density;  
  double K;        
  double G;        
};
 
template <typename T>
auto overstress(double viscosity, T accumulated_plastic_strain_rate)
{
  return viscosity * accumulated_plastic_strain_rate;
}
 
struct LinearHardening {
  double sigma_y;  
  double Hi;       
  double eta;      
 
  template <typename T1, typename T2>
  auto operator()(T1 accumulated_plastic_strain, T2 accumulated_plastic_strain_rate) const
  {
    return sigma_y + Hi * accumulated_plastic_strain + overstress(eta, accumulated_plastic_strain_rate);
  };
};
 
struct PowerLawHardening {
  double sigma_y;  
  double n;        
  double eps0;     
  double eta;      
 
  template <typename T1, typename T2>
  auto operator()(T1 accumulated_plastic_strain, T2 accumulated_plastic_strain_rate) const
  {
    using std::pow;
    return sigma_y * pow(1.0 + accumulated_plastic_strain / eps0, 1.0 / n) +
           overstress(eta, accumulated_plastic_strain_rate);
  };
};
 
struct VoceHardening {
  double sigma_y;          
  double sigma_sat;        
  double strain_constant;  
  double eta;              
 
  template <typename T1, typename T2>
  auto operator()(T1 accumulated_plastic_strain, T2 accumulated_plastic_strain_rate) const
  {
    using std::exp;
    auto& epsilon_p = accumulated_plastic_strain;
    auto Y_eq = sigma_sat - (sigma_sat - sigma_y) * exp(-epsilon_p / strain_constant);
    auto& epsilon_p_dot = accumulated_plastic_strain_rate;
    auto Y_vis = overstress(eta, epsilon_p_dot);
    return Y_eq + Y_vis;
  };
};
 
template <typename HardeningType>
struct J2SmallStrain {
  static constexpr int dim = 3;         
  static constexpr double tol = 1e-10;  
 
  double E;                 
  double nu;                
  HardeningType hardening;  
  double Hk;                
  double density;           
 
  struct State {
    tensor<double, dim, dim> plastic_strain;  
    double accumulated_plastic_strain;        
  };
 
  template <typename T>
  auto operator()(State& state, double dt, const tensor<T, dim, dim>& 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 sigma_b = 2.0 / 3.0 * Hk * state.plastic_strain;
    auto eta = s - sigma_b;
    auto q = sqrt(1.5) * norm(eta);
 
    // (ii) admissibility
    const double eqps_old = state.accumulated_plastic_strain;
    auto residual = [eqps_old, G, dt, *this](auto delta_eqps, auto trial_q) {
      auto eqps_dot = dt > 0.0 ? delta_eqps / dt : 0.0 * delta_eqps;
      return trial_q - (3.0 * G + Hk) * delta_eqps - this->hardening(eqps_old + delta_eqps, eqps_dot);
    };
    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, 0.0)) / (3.0 * G + Hk);
      auto [delta_eqps, status] = solve_scalar_equation(residual, 0.0, lower_bound, upper_bound, opts, q);
 
      auto Np = 1.5 * eta / 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;
  }
};
 
template <typename HardeningType>
struct J2 {
  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> Fpinv = DenseIdentity<3>();  
    double accumulated_plastic_strain;                    
  };
 
  template <typename T>
  auto operator()(State& state, double dt, const tensor<T, dim, dim>& 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 F = du_dX + I;
    auto Fe = dot(F, state.Fpinv);
    auto Ee = 0.5 * log_symm(dot(transpose(Fe), Fe));
    // From this point until the state variable update, the algorithm exactly coincides with the
    // small strain one.
    auto p = K * tr(Ee);
    auto s = 2.0 * G * dev(Ee);
    double q_val = sqrt(1.5) * norm(get_value(s));
 
    // (ii) admissibility
    const double eqps_old = state.accumulated_plastic_strain;
    auto residual = [eqps_old, G, dt, *this](auto delta_eqps, auto trial_mises) {
      auto eqps_dot = dt > 0.0 ? delta_eqps / dt : 0.0 * delta_eqps;
      return trial_mises - 3.0 * G * delta_eqps - this->hardening(eqps_old + delta_eqps, eqps_dot);
    };
    if (residual(0.0, q_val) > 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 = (q_val - hardening(eqps_old, 0.0)) / (3.0 * G);
      // safe to compute derivative of mises stress now that we're in yielding branch
      auto q = sqrt(1.5) * norm(s);
      auto [delta_eqps, status] = solve_scalar_equation(residual, 0.0, lower_bound, upper_bound, opts, q);
 
      if (!status.converged) {
        SLIC_WARNING("J2 solve failed to converge.");
      }
 
      auto Np = 1.5 * s / q;
 
      s = s - 2.0 * G * delta_eqps * Np;
      auto A = exp_symm(-delta_eqps * Np);
 
      auto Fpinv = dot(state.Fpinv, A);
 
      state.accumulated_plastic_strain += get_value(delta_eqps);
      state.Fpinv = get_value(Fpinv);
      // Mandel stress
      auto M = s + p * I;
      // convert to Piola
      Fe = dot(Fe, A);
      auto P = transpose(dot(dot(Fpinv, M), inv(Fe)));
      return P;
    } else {
      // I want to unify this branch with the yielding branch, but I'd need to declare Fpinv above,
      // and I don't know how to declare its type. BT 11/5/24
      // Mandel stress
      auto M = s + p * I;
      // convert to Piola
      auto P = transpose(dot(dot(state.Fpinv, M), inv(Fe)));
      return P;
    }
  }
};
 
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_;
  }
};
 
struct NeoHookeanWithFieldDensity {
  using State = Empty;  
 
  template <typename T, int dim, typename Density>
  SERAC_HOST_DEVICE auto pkStress(State& /* state */, const tensor<T, dim, dim>& du_dX, const Density&) const
  {
    using std::log1p;
    constexpr auto I = Identity<dim>();
    auto lambda = K - (2.0 / 3.0) * G;
    auto B_minus_I = dot(du_dX, transpose(du_dX)) + transpose(du_dX) + du_dX;
 
    auto logJ = log1p(detApIm1(du_dX));
    // Kirchoff stress, in form that avoids cancellation error when F is near I
    auto TK = lambda * logJ * I + G * B_minus_I;
 
    // Pull back to Piola
    auto F = du_dX + I;
    return dot(TK, inv(transpose(F)));
  }
 
  template <typename Density>
  SERAC_HOST_DEVICE auto density(const Density& density) const
  {
    return get<VALUE>(density);
  }
 
  double K;  
  double G;  
};
 
}  // namespace serac::solid_mechanics