green_saint_venant_thermoelastic.hpp

// 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)
 
#include "serac/numerics/functional/tensor.hpp"
#include "serac/numerics/functional/tuple.hpp"
 
namespace serac {
 
template <typename T>
auto greenStrain(const tensor<T, 3, 3>& 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 k;          
 
  struct State {
    double strain_trace;  
  };
 
  template <typename T1, typename T2, typename T3>
  auto operator()(State& state, const tensor<T1, 3, 3>& grad_u, T2 theta, const tensor<T3, 3>& 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<3>();
    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 P     = dot(F, S);
    const auto sigma = dot(P, transpose(F)) / det(F);
 
    // internal heat source
    const auto s0 = -3 * K * alpha * theta * (trEg - state.strain_trace);
 
    // heat flux
    const auto q0 = -k * grad_theta;
 
    state.strain_trace = get_value(trEg);
 
    return serac::tuple{sigma, C_v, s0, q0};
  }
 
  template <typename T1, typename T2>
  auto calculateFreeEnergy(const tensor<T1, 3, 3>& 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 k;          
 
  struct State {
    double strain_trace;  
  };
 
  template <typename T1, typename T2, typename T3, typename T4>
  auto operator()(State& state, const tensor<T1, 3, 3>& grad_u, T2 theta, const tensor<T3, 3>& 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<3>();
    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 P     = dot(F, S);
    const auto sigma = (dot(P, transpose(F))) / det(F);
 
    // internal heat source
    const auto s0 = -3 * K * alpha * theta * (trEg - state.strain_trace);
 
    // heat flux
    const auto q0 = -k * grad_theta;
 
    state.strain_trace = get_value(trEg);
 
    return serac::tuple{sigma, C_v, s0, q0};
  }
 
  template <typename T1, typename T2, typename T3>
  auto calculateFreeEnergy(const tensor<T1, 3, 3>& 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 serac