creepdam.cpp

Finite element program for mechanical problem. It can solve various problem in solid problem

#include "creepdam.h"
#include "matrix.h"
#include "vector.h"
#include "elastisomat.h"
#include "global.h"
#include "intpoints.h"
#include "vecttens.h"
#include <math.h>

/**
  This constructor inializes attributes to zero values.
*/
creepdam::creepdam (void)
{

}

/**
  This destructor is only for the formal purposes.
*/
creepdam::~creepdam (void)
{

}

/**
  This function computes material stiffnes %matrix.

  @param d - allocated matrix structure for material stiffness %matrix
  @param ipp - integration point number
  @param im - material index
  @param ido - index of internal variables for given material in the ipp eqother array


*/
void creepdam::matstiff (matrix &d,long ipp,long im, long ido)
{
  long ncompo = Mm->givencompeqother(ipp, im+1); 
  Mm->matstiff(d, ipp, im+2, ido+ncompo);
}


/**
  This function computes correct stresses in the integration point and stores
  them into ip stress array. It was used for Hinkley computations and supposes creep 
  strains to be read from backup file.
  
  @param ipp - integration point pointer
  7.10.2001

void creepdam::nlstresses (long ipp)
{
  long i;
  long ncomp = Mm->ip[ipp].ncompstr;
  vector epsback(ncomp);


  //  backup of the total strains
  for (i=0;i<ncomp;i++){
    epsback[i] = Mm->ip[ipp].strain[i];
  }

  // compute creep - actually it is supposed the creep strains are read from backup file
  //                 to the eqother array.
  //  initial values of elastic strain for damage = total strain - creep strain
  for (i=0;i<ncomp;i++){
    Mm->ip[ipp].strain[i] -= Mm->ip[ipp].eqother[i];
  }
  
  // compute damage, it has index 1, eqother values for damage starts at ncomp index
  Mm->computenlstresses(ipp, 1, 2*ncomp);
  
  //  recovery of the total strains
  for (i=0;i<ncomp;i++){
    Mm->ip[ipp].strain[i] = epsback[i];
  }

}
*/

/**
  This function computes correct stresses in the integration point and stores
  them into ip stress array.
  
  @param ipp - integration point pointer
  @param im  - index of the material in the ipp tm and idm arrays. The standard value is zero.
  7.10.2001
*/
void creepdam::nlstresses (long ipp, long im, long ido)
{
  long i, idem, nm;
  long ncomp = Mm->ip[ipp].ncompstr;
  long ncompo;
  double nu;
  vector epsback(ncomp);
  vector sigback(ncomp);
  vector epscr(ncomp);


  if (Mm->ip[ipp].ssst == planestress)
  {
    idem = Mm->ip[ipp].gemid();
    nu = Mm->eliso[Mm->ip[ipp].idm[idem]].nu;
    Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (Mm->ip[ipp].strain[0]+Mm->ip[ipp].strain[1]);
  }

  if (Mp->phase == 1)
  {
    // compute creep , it has index1, creep eqother values starts from 0.
    Mm->computenlstresses(ipp, 1, ido+0);
    //  backup of stress increments from creep
    for (i=0;i<ncomp;i++)
      sigback[i] = Mm->ip[ipp].stress[i];
    // compute increment of stresses due to temperature strain
    nm=Mm->ip[ipp].nm-1;
    ncompo = Mm->givencompeqother(ipp, 0);
    ncompo -= Mm->givencompeqother(ipp, nm);
    Mm->computenlstresses(ipp, nm, ncompo);
    //  compute total stress increment
    for (i=0;i<ncomp;i++)
      Mm->ip[ipp].stress[i] += sigback[i];
  }  
  
  if (Mp->phase == 2)
  {
    //  backup of the total strains
    for (i=0;i<ncomp;i++)
      epsback[i] = Mm->ip[ipp].strain[i];

    //  initial values of elastic strain for damage = total strain - creep strain
    Mm->giveirrstrains(ipp, im+1, ido+0, epscr);
    for (i=0;i<ncomp;i++)
      Mm->ip[ipp].strain[i] -= epscr[i];

    // number of components in eqother array for creep 
    ncompo = Mm->givencompeqother(ipp, im+1); 
    // compute damage, it has index 2, eqother values for damage starts at ncompo index
    Mm->computenlstresses(ipp, im+2, ido+ncompo);
  
    //  recovery of the total strains and total stress increment
    for (i=0;i<ncomp;i++)
      Mm->ip[ipp].strain[i] = epsback[i];

    // actualization of temperature model
    if ((im == 0) && (Mm->ip[ipp].hmt & 1))
    {
      nm=Mm->ip[ipp].nm-1;
      ncompo = Mm->givencompeqother(ipp, 0);
      ncompo -= Mm->givencompeqother(ipp, nm);
      // actualization of previous temperature
      Mm->computenlstresses(ipp, nm, ncompo);
    }
  }
}


/**
  This function updates values in the other array reached in the previous equlibrium state to
  values reached in the new actual equilibrium state.

  @param ipp - integration point number in the mechmat ip array.
  @param im - material index
  @param ido - index of internal variables for given material in the ipp other array

*/
void creepdam::updateval (long ipp, long im, long ido)
{
  //  update values of plasticity model
  //  long ncompo = 2*Mm->ip[ipp].ncompstr; // for Hinkley
  //  Mm->updateipvalmat (ipp,1,ncompo); // for Hinkley
  long ncompo = Mm->givencompeqother(ipp, im+1); 
  Mm->updateipvalmat (ipp,1,ido+0); 
  Mm->updateipvalmat (ipp,2,ido+ncompo);  
}



/**
   Function initializes eqother array with initial initial values.
   Actual value of array Mm->tempr[ipp] is taken as initial temperature.
   
   @param ipp - integration point pointer
   @param im - material index
   @param ido - index of internal variables for given material in the ipp other array

   7.6.2005, TKo
*/
void creepdam::initvalues (long ipp, long im, long ido)
{
  long ncompo = Mm->givencompeqother(ipp, im+1); 
  Mm->initvalues(ipp,im+1,0); 
  Mm->initvalues(ipp,im+2,ncompo); 
}


/**
  This function returns the value of tensile strength

  @param ipp - integration point number in the mechmat ip array.
  @param im - material index
  @param ido - index of internal variables for given material in the ipp other array
  
  TKo
*/
double creepdam::give_actual_ft (long ipp, long im, long ido)
{
  double ft;
  ft = Mm->give_actual_ft(ipp, im+1, ido);
  return ft;
}