green_saint_venant_thermoelastic.hpp

// Copyright (c) Lawrence Livermore National Security, LLC and
// other Smith Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include "smith/numerics/functional/tensor.hpp"
#include "smith/numerics/functional/tuple.hpp"
 
namespace smith::thermomechanics {
 
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));
}
 
template <typename EType>
auto bulkModulus(EType E, double nu)
{
  return E / (3.0 * (1.0 - 2.0 * nu));
}
 
template <typename EType>
auto shearModulus(EType E, double nu)
{
  return 0.5 * E / (1.0 + nu);
}
 
template <typename EType, typename AlphaType, typename TGradU, typename TTheta, int dim>
auto greenSaintVenantPiola(EType E, double nu, AlphaType alpha, double theta_ref,
                           const tensor<TGradU, dim, dim>& grad_u, TTheta theta)
{
  const auto K = bulkModulus(E, nu);
  const auto G = shearModulus(E, nu);
  static constexpr auto I = Identity<dim>();
  auto F = grad_u + I;
  const auto Eg = greenStrain(grad_u);
  const auto trEg = tr(Eg);
  const auto S = 2.0 * G * dev(Eg) + K * (trEg - static_cast<double>(dim) * alpha * (theta - theta_ref)) * I;
  return dot(F, S);
}
 
template <typename TGradTheta, int dim>
auto fourierHeatFlux(double kappa, const tensor<TGradTheta, dim>& grad_theta)
{
  return -kappa * grad_theta;
}
 
struct GreenSaintVenantThermoelasticMaterial {
  double density;    
  double E;          
  double nu;         
  double C_v;        
  double alpha;      
  double theta_ref;  
  double kappa;      
 
  struct State {
    double strain_trace;  
  };
 
  template <typename T1, typename T2, typename T3, int dim>
  auto operator()(State& state, const tensor<T1, dim, dim>& grad_u, T2 theta, const tensor<T3, dim>& grad_theta) const
  {
    const auto Eg = greenStrain(grad_u);
    const auto trEg = tr(Eg);
    const auto Piola = greenSaintVenantPiola(E, nu, alpha, theta_ref, grad_u, theta);
    const auto s0 = -3.0 * bulkModulus(E, nu) * alpha * theta * (trEg - state.strain_trace);
    const auto q0 = fourierHeatFlux(kappa, grad_theta);
 
    state.strain_trace = get_value(trEg);
 
    return smith::tuple{Piola, C_v, s0, q0};
  }
 
  template <typename T1, typename T2, int dim>
  auto calculateFreeEnergy(const tensor<T1, dim, dim>& grad_u, T2 theta) const
  {
    const double K = E / (3.0 * (1.0 - 2.0 * nu));
    const double G = 0.5 * E / (1.0 + nu);
    auto strain = greenStrain(grad_u);
    auto trE = tr(strain);
    auto psi_1 = G * squared_norm(dev(strain)) + 0.5 * K * trE * trE;
    using std::log;
    auto logT = log(theta / theta_ref);
    auto psi_2 = C_v * (theta - theta_ref - theta * logT);
    auto psi_3 = -3.0 * K * alpha * (theta - theta_ref) * trE;
    return psi_1 + psi_2 + psi_3;
  }
};
 
struct ParameterizedGreenSaintVenantThermoelasticMaterial {
  double density;    
  double E;          
  double nu;         
  double C_v;        
  double alpha0;     
  double theta_ref;  
  double kappa;      
 
  struct State {
    double strain_trace;  
  };
 
  template <typename T1, typename T2, typename T3, typename T4, int dim>
  auto operator()(State& state, const tensor<T1, dim, dim>& grad_u, T2 theta, const tensor<T3, dim>& grad_theta,
                  T4 thermal_expansion_scaling) const
  {
    auto [scale, unused] = thermal_expansion_scaling;
    const auto Eg = greenStrain(grad_u);
    const auto trEg = tr(Eg);
    auto alpha = alpha0 * scale;
    const auto Piola = greenSaintVenantPiola(E, nu, alpha, theta_ref, grad_u, theta);
    const auto s0 = -3.0 * bulkModulus(E, nu) * alpha * theta * (trEg - state.strain_trace);
    const auto q0 = fourierHeatFlux(kappa, grad_theta);
 
    state.strain_trace = get_value(trEg);
 
    return smith::tuple{Piola, C_v, s0, q0};
  }
 
  template <typename T1, typename T2, typename T3, int dim>
  auto calculateFreeEnergy(const tensor<T1, dim, dim>& grad_u, T2 theta, T3 thermal_expansion_scaling) const
  {
    auto [scale, unused] = thermal_expansion_scaling;
    const double K = bulkModulus(E, nu);
    const double G = shearModulus(E, nu);
    auto strain = greenStrain(grad_u);
    auto trE = tr(strain);
    const auto alpha = alpha0 * scale;
    auto psi_1 = G * squared_norm(dev(strain)) + 0.5 * K * trE * trE;
    using std::log;
    auto logT = log(theta / theta_ref);
    auto psi_2 = C_v * (theta - theta_ref - theta * logT);
    auto psi_3 = -3.0 * K * alpha * (theta - theta_ref) * trE;
    return psi_1 + psi_2 + psi_3;
  }
};
 
}  // namespace smith::thermomechanics