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));
}
 
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 double K = E / (3.0 * (1.0 - 2.0 * nu));
    const double G = 0.5 * E / (1.0 + nu);
    static constexpr auto I = Identity<dim>();
    auto F = grad_u + I;
    const auto Eg = greenStrain(grad_u);
    const auto trEg = tr(Eg);
 
    // stress
    const auto S = 2.0 * G * dev(Eg) + K * (trEg - 3.0 * alpha * (theta - theta_ref)) * I;
    const auto Piola = dot(F, S);
 
    // internal heat source
    const auto s0 = -3.0 * K * alpha * theta * (trEg - state.strain_trace);
 
    // heat flux
    const auto q0 = -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 double K = E / (3.0 * (1.0 - 2.0 * nu));
    const double G = 0.5 * E / (1.0 + nu);
    static constexpr auto I = Identity<dim>();
    auto F = grad_u + I;
    const auto Eg = greenStrain(grad_u);
    const auto trEg = tr(Eg);
    auto alpha = alpha0 * scale;
 
    // stress
    const auto S = 2.0 * G * dev(Eg) + K * (trEg - 3.0 * alpha * (theta - theta_ref)) * I;
    const auto Piola = dot(F, S);
 
    // internal heat source
    const auto s0 = -3.0 * K * alpha * theta * (trEg - state.strain_trace);
 
    // heat flux
    const auto q0 = -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 = 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);
    const double 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