oomph::IsotropicStrainEnergyFunctionConstitutiveLaw Class Reference

#include <constitutive_laws.h>

Inheritance diagram for oomph::IsotropicStrainEnergyFunctionConstitutiveLaw:

Public Member Functions
	IsotropicStrainEnergyFunctionConstitutiveLaw (StrainEnergyFunction *const &strain_energy_function_pt)
	Constructor takes a pointer to the strain energy function. More...
void	calculate_second_piola_kirchhoff_stress (const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
void	calculate_second_piola_kirchhoff_stress (const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma_dev, DenseMatrix< double > &G_contra, double &Gdet)
void	calculate_second_piola_kirchhoff_stress (const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma_dev, DenseMatrix< double > &Gcontra, double &gen_dil, double &inv_kappa)
bool	requires_incompressibility_constraint ()
Public Member Functions inherited from oomph::ConstitutiveLaw
	ConstitutiveLaw ()
	Empty constructor. More...
virtual	~ConstitutiveLaw ()
	Empty virtual destructor. More...
virtual void	calculate_d_second_piola_kirchhoff_stress_dG (const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, RankFourTensor< double > &d_sigma_dG, const bool &symmetrize_tensor=true)
virtual void	calculate_d_second_piola_kirchhoff_stress_dG (const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, const double &detG, const double &interpolated_solid_p, RankFourTensor< double > &d_sigma_dG, DenseMatrix< double > &d_detG_dG, const bool &symmetrize_tensor=true)
virtual void	calculate_d_second_piola_kirchhoff_stress_dG (const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, const double &gen_dil, const double &inv_kappa, const double &interpolated_solid_p, RankFourTensor< double > &d_sigma_dG, DenseMatrix< double > &d_gen_dil_dG, const bool &symmetrize_tensor=true)

Private Attributes
StrainEnergyFunction *	Strain_energy_function_pt
	Pointer to the strain energy function. More...

Additional Inherited Members
Protected Member Functions inherited from oomph::ConstitutiveLaw
bool	is_matrix_square (const DenseMatrix< double > &M)
	Test whether a matrix is square. More...
bool	are_matrices_of_equal_dimensions (const DenseMatrix< double > &M1, const DenseMatrix< double > &M2)
	Test whether two matrices are of equal dimensions. More...
void	error_checking_in_input (const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
double	calculate_contravariant (const DenseMatrix< double > &Gcov, DenseMatrix< double > &Gcontra)
	The function to calculate the contravariant tensor from a covariant one. More...
void	calculate_d_contravariant_dG (const DenseMatrix< double > &Gcov, RankFourTensor< double > &dGcontra_dG, DenseMatrix< double > &d_detG_dG)

Detailed Description

////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// A class for constitutive laws derived from strain-energy functions. Theory is in Green and Zerna.

Constructor & Destructor Documentation

IsotropicStrainEnergyFunctionConstitutiveLaw()

oomph::IsotropicStrainEnergyFunctionConstitutiveLaw::IsotropicStrainEnergyFunctionConstitutiveLaw ( StrainEnergyFunction *const &  strain_energy_function_pt )

inline

Constructor takes a pointer to the strain energy function.

      : ConstitutiveLaw(), Strain_energy_function_pt(strain_energy_function_pt)
    {
    }

Member Function Documentation

calculate_second_piola_kirchhoff_stress() [1/3]

void oomph::IsotropicStrainEnergyFunctionConstitutiveLaw::calculate_second_piola_kirchhoff_stress ( const DenseMatrix< double > &  g, const DenseMatrix< double > &  G, DenseMatrix< double > &  sigma  )

virtual

Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed and deformed metric tensor and the matrix in which to return the stress tensor. Uses correct 3D invariants for 2D (plane strain) problems.

////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed and deformed metric tensor and the matrix in which to return the stress tensor. Uses correct 3D invariants for 2D (plane strain) problems.

Implements oomph::ConstitutiveLaw.

  {
// Error checking
#ifdef PARANOID
    error_checking_in_input(g, G, sigma);
#endif
 
    // Find the dimension of the problem
    unsigned dim = g.nrow();
 
#ifdef PARANOID
    if (dim == 1)
    {
      std::string function_name =
        "IsotropicStrainEnergyFunctionConstitutiveLaw::";
      function_name += "calculate_second_piola_kirchhoff_stress()";
 
      throw OomphLibError("Check constitutive equations carefully when dim=1",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Calculate the contravariant undeformed and deformed metric tensors
    // and get the determinants of the metric tensors
    DenseMatrix<double> gup(dim), Gup(dim);
    double detg = calculate_contravariant(g, gup);
    double detG = calculate_contravariant(G, Gup);
 
    // Calculate the strain invariants
    Vector<double> I(3, 0.0);
    // The third strain invaraint is the volumetric change
    I[2] = detG / detg;
    // The first and second are a bit more complex --- see G&Z
    for (unsigned i = 0; i < dim; i++)
    {
      for (unsigned j = 0; j < dim; j++)
      {
        I[0] += gup(i, j) * G(i, j);
        I[1] += g(i, j) * Gup(i, j);
      }
    }
 
    // If 2D we assume plane strain: In this case the 3D tensors have
    // a 1 on the diagonal and zeroes in the off-diagonals of their
    // third rows and columns. Only effect: Increase the first two
    // invariants by one; rest of the computation can just be performed
    // over the 2d set of coordinates.
    if (dim == 2)
    {
      I[0] += 1.0;
      I[1] += 1.0;
    }
 
    // Second strain invariant is multiplied by the third.
    I[1] *= I[2];
 
    // Calculate the derivatives of the strain energy function wrt the
    // strain invariants
    Vector<double> dWdI(3, 0.0);
    Strain_energy_function_pt->derivatives(I, dWdI);
 
 
    // Only bother to compute the tensor B^{ij} (Green & Zerna notation)
    // if the derivative wrt the second strain invariant is non-zero
    DenseMatrix<double> Bup(dim, dim, 0.0);
    if (std::fabs(dWdI[1]) > 0.0)
    {
      for (unsigned i = 0; i < dim; i++)
      {
        for (unsigned j = 0; j < dim; j++)
        {
          Bup(i, j) = I[0] * gup(i, j);
          for (unsigned r = 0; r < dim; r++)
          {
            for (unsigned s = 0; s < dim; s++)
            {
              Bup(i, j) -= gup(i, r) * gup(j, s) * G(r, s);
            }
          }
        }
      }
    }
 
    // Now set the values of the functions phi, psi and p (Green & Zerna
    // notation) Note that the Green & Zerna stress \tau^{ij} is
    // s^{ij}/sqrt(I[2]), where s^{ij} is the desired second Piola-Kirchhoff
    // stress tensor so we multiply their constants by sqrt(I[2])
    double phi = 2.0 * dWdI[0];
    double psi = 2.0 * dWdI[1];
    double p = 2.0 * dWdI[2] * I[2];
 
    // Put it all together to get the stress
    for (unsigned i = 0; i < dim; i++)
    {
      for (unsigned j = 0; j < dim; j++)
      {
        sigma(i, j) = phi * gup(i, j) + psi * Bup(i, j) + p * Gup(i, j);
      }
    }
  }

References oomph::ConstitutiveLaw::calculate_contravariant(), oomph::StrainEnergyFunction::derivatives(), oomph::ConstitutiveLaw::error_checking_in_input(), boost::multiprecision::fabs(), G, i, I, j, oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p, UniformPSDSelfTest::r, s, calibrate::sigma, Strain_energy_function_pt, and oomph::Global_string_for_annotation::string().

calculate_second_piola_kirchhoff_stress() [2/3]

void oomph::IsotropicStrainEnergyFunctionConstitutiveLaw::calculate_second_piola_kirchhoff_stress ( const DenseMatrix< double > &  g, const DenseMatrix< double > &  G, DenseMatrix< double > &  sigma_dev, DenseMatrix< double > &  Gup, double &  detG  )

virtual

Calculate the deviatoric part \( \overline{ \sigma^{ij}}\) of the contravariant 2nd Piola Kirchhoff stress tensor \( \sigma^{ij}\). Also return the contravariant deformed metric tensor and the determinant of the deformed metric tensor. This form is appropriate for truly-incompressible materials for which \( \sigma^{ij} = - p G^{ij} +\overline{ \sigma^{ij}} \) where the "pressure" \( p \) is determined by \( \det G_{ij} - \det g_{ij} = 0 \).

Calculate the deviatoric part \( \overline{ \sigma^{ij}}\) of the contravariant 2nd Piola Kirchhoff stress tensor \( \sigma^{ij}\). Also return the contravariant deformed metric tensor and the determinant of the deformed metric tensor. Uses correct 3D invariants for 2D (plane strain) problems. This is the version for the pure incompressible formulation.

Reimplemented from oomph::ConstitutiveLaw.

  {
// Error checking
#ifdef PARANOID
    error_checking_in_input(g, G, sigma_dev);
#endif
 
    // Find the dimension of the problem
    unsigned dim = g.nrow();
 
 
#ifdef PARANOID
    if (dim == 1)
    {
      std::string function_name =
        "IsotropicStrainEnergyFunctionConstitutiveLaw::";
      function_name += "calculate_second_piola_kirchhoff_stress()";
 
      throw OomphLibError("Check constitutive equations carefully when dim=1",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Calculate the contravariant undeformed and deformed metric tensors
    DenseMatrix<double> gup(dim);
    // Don't need this determinant
    (void)calculate_contravariant(g, gup);
    // These are passed back
    detG = calculate_contravariant(G, Gup);
 
    // Calculate the strain invariants
    Vector<double> I(3, 0.0);
    // The third strain invaraint must be one (incompressibility)
    I[2] = 1.0;
    // The first and second are a bit more complex
    for (unsigned i = 0; i < dim; i++)
    {
      for (unsigned j = 0; j < dim; j++)
      {
        I[0] += gup(i, j) * G(i, j);
        I[1] += g(i, j) * Gup(i, j);
      }
    }
 
    // If 2D we assume plane strain: In this case the 3D tensors have
    // a 1 on the diagonal and zeroes in the off-diagonals of their
    // third rows and columns. Only effect: Increase the first two
    // invariants by one; rest of the computation can just be performed
    // over the 2d set of coordinates.
    if (dim == 2)
    {
      I[0] += 1.0;
      I[1] += 1.0;
    }
 
    // Calculate the derivatives of the strain energy function wrt the
    // strain invariants
    Vector<double> dWdI(3, 0.0);
    Strain_energy_function_pt->derivatives(I, dWdI);
 
    // Only bother to compute the tensor B^{ij} (Green & Zerna notation)
    // if the derivative wrt the second strain invariant is non-zero
    DenseMatrix<double> Bup(dim, dim, 0.0);
    if (std::fabs(dWdI[1]) > 0.0)
    {
      for (unsigned i = 0; i < dim; i++)
      {
        for (unsigned j = 0; j < dim; j++)
        {
          Bup(i, j) = I[0] * gup(i, j);
          for (unsigned r = 0; r < dim; r++)
          {
            for (unsigned s = 0; s < dim; s++)
            {
              Bup(i, j) -= gup(i, r) * gup(j, s) * G(r, s);
            }
          }
        }
      }
    }
 
    // Now set the values of the functions phi and psi (Green & Zerna notation)
    double phi = 2.0 * dWdI[0];
    double psi = 2.0 * dWdI[1];
    // Calculate the trace/dim of the first two terms of the stress tensor
    // phi g^{ij} + psi B^{ij} (see Green & Zerna)
    double K;
    // In two-d, we cannot use the strain invariants directly
    // but can use symmetry of the tensors involved
    if (dim == 2)
    {
      K = 0.5 * ((I[0] - 1.0) * phi +
                 psi * (Bup(0, 0) * G(0, 0) + Bup(1, 1) * G(1, 1) +
                        2.0 * Bup(0, 1) * G(0, 1)));
    }
    // In three-d we can make use of the strain invariants, see Green & Zerna
    else
    {
      K = (I[0] * phi + 2.0 * I[1] * psi) / 3.0;
    }
 
    // Put it all together to get the stress, subtracting the trace of the
    // first two terms to ensure that the stress is deviatoric, which means
    // that the computed pressure is the mechanical pressure
    for (unsigned i = 0; i < dim; i++)
    {
      for (unsigned j = 0; j < dim; j++)
      {
        sigma_dev(i, j) = phi * gup(i, j) + psi * Bup(i, j) - K * Gup(i, j);
      }
    }
  }

References oomph::ConstitutiveLaw::calculate_contravariant(), oomph::StrainEnergyFunction::derivatives(), oomph::ConstitutiveLaw::error_checking_in_input(), boost::multiprecision::fabs(), G, i, I, j, PlanarWave::K, oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, UniformPSDSelfTest::r, s, Strain_energy_function_pt, and oomph::Global_string_for_annotation::string().

calculate_second_piola_kirchhoff_stress() [3/3]

void oomph::IsotropicStrainEnergyFunctionConstitutiveLaw::calculate_second_piola_kirchhoff_stress ( const DenseMatrix< double > &  g, const DenseMatrix< double > &  G, DenseMatrix< double > &  sigma_dev, DenseMatrix< double > &  Gup, double &  gen_dil, double &  inv_kappa  )

virtual

Calculate the deviatoric part of the contravariant 2nd Piola Kirchoff stress tensor. Also return the contravariant deformed metric tensor, the generalised dilatation, \( d, \) and the inverse of the bulk modulus \( \kappa\). This form is appropriate for near-incompressible materials for which \( \sigma^{ij} = -p G^{ij} + \overline{ \sigma^{ij}} \) where the "pressure" \( p \) is determined from \( p / \kappa - d =0 \).

Calculate the deviatoric part of the contravariant 2nd Piola Kirchoff stress tensor. Also return the contravariant deformed metric tensor, the generalised dilatation, \( d, \) and the inverse of the bulk modulus \( \kappa\). Uses correct 3D invariants for 2D (plane strain) problems. This is the version for the near-incompressible formulation.

Reimplemented from oomph::ConstitutiveLaw.

  {
// Error checking
#ifdef PARANOID
    error_checking_in_input(g, G, sigma_dev);
#endif
 
    // Find the dimension of the problem
    unsigned dim = g.nrow();
 
#ifdef PARANOID
    if (dim == 1)
    {
      std::string function_name =
        "IsotropicStrainEnergyFunctionConstitutiveLaw::";
      function_name += "calculate_second_piola_kirchhoff_stress()";
 
      throw OomphLibError("Check constitutive equations carefully when dim=1",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Calculate the contravariant undeformed and deformed metric tensors
    // and get the determinants of the metric tensors
    DenseMatrix<double> gup(dim);
    double detg = calculate_contravariant(g, gup);
    double detG = calculate_contravariant(G, Gup);
 
    // Calculate the strain invariants
    Vector<double> I(3, 0.0);
    // The third strain invaraint is the volumetric change
    I[2] = detG / detg;
    // The first and second are a bit more complex --- see G&Z
    for (unsigned i = 0; i < dim; i++)
    {
      for (unsigned j = 0; j < dim; j++)
      {
        I[0] += gup(i, j) * G(i, j);
        I[1] += g(i, j) * Gup(i, j);
      }
    }
 
    // If 2D we assume plane strain: In this case the 3D tensors have
    // a 1 on the diagonal and zeroes in the off-diagonals of their
    // third rows and columns. Only effect: Increase the first two
    // invariants by one; rest of the computation can just be performed
    // over the 2d set of coordinates.
    if (dim == 2)
    {
      I[0] += 1.0;
      I[1] += 1.0;
    }
 
    // Second strain invariant is multiplied by the third.
    I[1] *= I[2];
 
    // Calculate the derivatives of the strain energy function wrt the
    // strain invariants
    Vector<double> dWdI(3, 0.0);
    Strain_energy_function_pt->derivatives(I, dWdI);
 
    // Only bother to calculate the tensor B^{ij} (Green & Zerna notation)
    // if the derivative wrt the second strain invariant is non-zero
    DenseMatrix<double> Bup(dim, dim, 0.0);
    if (std::fabs(dWdI[1]) > 0.0)
    {
      for (unsigned i = 0; i < dim; i++)
      {
        for (unsigned j = 0; j < dim; j++)
        {
          Bup(i, j) = I[0] * gup(i, j);
          for (unsigned r = 0; r < dim; r++)
          {
            for (unsigned s = 0; s < dim; s++)
            {
              Bup(i, j) -= gup(i, r) * gup(j, s) * G(r, s);
            }
          }
        }
      }
    }
 
    // Now set the values of the functions phi and psi (Green & Zerna notation)
    // but multiplied by sqrt(I[2]) to recover the second Piola-Kirchhoff stress
    double phi = 2.0 * dWdI[0];
    double psi = 2.0 * dWdI[1];
 
    // Calculate the trace/dim of the first two terms of the stress tensor
    // phi g^{ij} + psi B^{ij} (see Green & Zerna)
    double K;
    // In two-d, we cannot use the strain invariants directly,
    // but we can use symmetry of the tensors involved
    if (dim == 2)
    {
      K = 0.5 * ((I[0] - 1.0) * phi +
                 psi * (Bup(0, 0) * G(0, 0) + Bup(1, 1) * G(1, 1) +
                        2.0 * Bup(0, 1) * G(0, 1)));
    }
    // In three-d we can make use of the strain invariants
    else
    {
      K = (I[0] * phi + 2.0 * I[1] * psi) / 3.0;
    }
 
    // Choose inverse kappa to be one...
    inv_kappa = 1.0;
 
    //...then the generalised dilation is the same as p  in Green & Zerna's
    // notation, but multiplied by sqrt(I[2]) with the addition of the
    // terms that are subtracted to make the other part of the stress
    // deviatoric
    gen_dil = 2.0 * dWdI[2] * I[2] + K;
 
    // Calculate the deviatoric part of the stress by subtracting
    // the computed trace/dim
    for (unsigned i = 0; i < dim; i++)
    {
      for (unsigned j = 0; j < dim; j++)
      {
        sigma_dev(i, j) = phi * gup(i, j) + psi * Bup(i, j) - K * Gup(i, j);
      }
    }
  }

References oomph::ConstitutiveLaw::calculate_contravariant(), oomph::StrainEnergyFunction::derivatives(), oomph::ConstitutiveLaw::error_checking_in_input(), boost::multiprecision::fabs(), G, i, I, j, PlanarWave::K, oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, UniformPSDSelfTest::r, s, Strain_energy_function_pt, and oomph::Global_string_for_annotation::string().

requires_incompressibility_constraint()

bool oomph::IsotropicStrainEnergyFunctionConstitutiveLaw::requires_incompressibility_constraint ( )

inlinevirtual

State if the constitutive equation requires an incompressible formulation in which the volume constraint is enforced explicitly. Used as a sanity check in PARANOID mode. This is determined by interrogating the associated strain energy function.

Implements oomph::ConstitutiveLaw.

    {
      return Strain_energy_function_pt->requires_incompressibility_constraint();
    }

References oomph::StrainEnergyFunction::requires_incompressibility_constraint(), and Strain_energy_function_pt.

Member Data Documentation

Strain_energy_function_pt

StrainEnergyFunction* oomph::IsotropicStrainEnergyFunctionConstitutiveLaw::Strain_energy_function_pt

private

Pointer to the strain energy function.

Referenced by calculate_second_piola_kirchhoff_stress(), and requires_incompressibility_constraint().

---

The documentation for this class was generated from the following files:

• constitutive_laws.h
• constitutive_laws.cc