AnalyticalSolutions.hpp

Go to the documentation of this file.

#ifndef __ANALYTICAL_SOLUTIONS_HPP__
#define __ANALYTICAL_SOLUTIONS_HPP__
 
#include <cmath>
 
namespace AnalyticalSolutions {
 
/**
 * \brief Lamé analytical solution for a hollow cylinder under radial pressure
 * with a linear isotropic Hooke material.
 *
 * The cylinder axis is aligned with z, while the solution is evaluated in the
 * x-y cross-section. The z coordinate is ignored. SPACE_DIM == 2 is the
 * plane-strain case with only the in-plane stress and displacement components
 * returned. For SPACE_DIM == 3 this uses the plane-strain axial stress and zero
 * axial displacement.
 */
template <int SPACE_DIM> struct HollowCylinderUnderRadialPressure {
 
  HollowCylinderUnderRadialPressure(
      const double inner_radius, const double outer_radius,
      const double inner_pressure, const double outer_pressure,
      const double young_modulus, const double poisson_ratio)
      : E(young_modulus), nu(poisson_ratio) {
    const double a2 = sqr(inner_radius);
    const double b2 = sqr(outer_radius);
    const double denom = b2 - a2;
 
    A = (inner_pressure * a2 - outer_pressure * b2) / denom;
    B = (inner_pressure - outer_pressure) * a2 * b2 / denom;
  }
 
  /**
   * \brief Return the Cauchy stress in symmetric tensor storage at point
   * (x, y, z).
   */
  MoFEM::MatrixDouble stress(const double x, const double y,
                             const double) const {
    MoFEM::MatrixDouble stress((SPACE_DIM * (SPACE_DIM + 1)) / 2, 1);
    stress.clear();
 
    const double r2 = sqr(x) + sqr(y);
    const double sigma_r = A - B / r2;
    const double sigma_t = A + B / r2;
    const double inv_r = 1. / std::sqrt(r2);
    const double c = x * inv_r;
    const double s = y * inv_r;
 
    auto t_stress = MoFEM::getFTensor2SymmetricFromMat<
        SPACE_DIM, -1,
        MoFEM::DataLayoutTraits<MoFEM::DataLayout::CoeffsByGauss>>(stress);
    t_stress(0, 0) = sigma_r * sqr(c) + sigma_t * sqr(s);
    t_stress(1, 1) = sigma_r * sqr(s) + sigma_t * sqr(c);
    t_stress(0, 1) = (sigma_r - sigma_t) * s * c;
    if constexpr (SPACE_DIM == 3)
      t_stress(2, 2) = 2. * nu * A;
 
    return stress;
  }
 
  /**
   * \brief Return the displacement vector at point (x, y, z).
   */
  MoFEM::VectorDouble displacement(const double x, const double y,
                                   const double) const {
    MoFEM::VectorDouble disp(SPACE_DIM);
    disp.clear();
 
    const double r2 = sqr(x) + sqr(y);
    const double r = std::sqrt(r2);
    const double u_r =
        (1. + nu) / E * ((1. - 2. * nu) * A * r + B / r);
 
    disp[0] = u_r * x / r;
    disp[1] = u_r * y / r;
    if constexpr (SPACE_DIM == 3)
      disp[2] = 0.;
 
    return disp;
  }
 
private:
  static constexpr double sqr(const double v) { return v * v; }
 
  const double E;
  const double nu;
  double A;
  double B;
};
 
} // namespace AnalyticalSolutions
 
#endif // __ANALYTICAL_SOLUTIONS_HPP__