anisodam.cpp

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

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

#define nijac 20
#define limit 1.0e-8
#ifndef M_PI
 #define M_PI 3.14159265358979323846
#endif





/**
  This constructor inializes attributes to zero values.
  
  TKo 9.2006
*/
anisodam::anisodam (void)
{
  cde = corr_off;
  ac = 0.0;  bc = 0.0;  y0c = 0.0;  at = 0.0;  bt = 0.0;  y0t = 0.0;
  a = 0.0;  b = 0.0;  y0 = 0.0; gfc = 0.0; gft = 0.0; gf = 0.0;
}

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

}

/**
  This function reads material parameters from the opened text file given
  by the parameter in.

  @param in - pointer to the opened input text file
  
  TKo 9.2006
*/
void anisodam::read (XFILE *in)
{
  xfscanf (in, "%d", (int *)&cde);
  switch (cde)
  {
    case corr_off:
      xfscanf (in, "%lf %lf %lf", &ac, &bc, &y0c);
      xfscanf (in, "%lf %lf %lf", &at, &bt, &y0t);
      xfscanf (in, "%lf %lf %lf", &a, &b, &y0);
      break;
    case corr_on:
      xfscanf (in, "%lf %lf %lf", &gfc , &bc, &y0c);
      xfscanf (in, "%lf %lf %lf", &gft, &bt, &y0t);
      xfscanf (in, "%lf %lf %lf", &gf, &b, &y0);
      break;
    default:
      fprintf(stderr, "\nUnknown type of correction of disipated energy is required\n");
      fprintf(stderr, "in function anisodam::read (file %s, line : %d)\n", __FILE__, __LINE__);
  }
}



/** 
    Function computes damage driving force for volumetric damage parameter.

    @param  ipp - integration point number
    @param  peps - %vector of principal strains
    
    @retval Function returns real value of damage driving force.

    TKo - 9.2006
*/
double anisodam::damdrvforce_vol(long ipp, vector &peps)
{
  long i;
  double ev = 0.0;
//  double e = Mm->give_actual_ym(ipp);
  long idem = Mm->ip[ipp].gemid();;
  long im = Mm->ip[ipp].idm[idem];
  double nu = Mm->eliso[im].nu;
  double k0;  
  // strain tensor due to positive principal strains
  for (i = 0; i < 3; i++)
    ev += peps(i);

  // take only positive volumetric strain  
  if (ev <= 0.0)
    return 0.0;

  // positive dimensionless volumetric damage driving force
  k0 = 1.0/(3.0*(1.0-2.0*nu));
  k0 *= 3.0/2.0*ev*ev;  
  return k0;
}

/** 
    Function computes damage driving force for deviatoric damage parameters.

    @param  ipp  - integration point number
    @param  peps - %vector of principal strains
    @param  t    - %matrix of principal directions stored in colomns for given
                   principal strains
    @param  pyt  - %vector of principal values of driving forces for tension damage parameters
    @param  pyc  - %vector of principal values of driving forces for compression damage parameters        

    TKo - 9.2006
*/
void anisodam::damdrvforce_dev(long ipp, vector &peps, matrix &t, vector &pyt, vector &pyc)
{
  long i;
  matrix et(3,3), ec(3,3), pe(3,3), tyt(3,3), tyc(3,3);  
//  double e = Mm->give_actual_ym(ipp);
  long idem = Mm->ip[ipp].gemid();;
  long im = Mm->ip[ipp].idm[idem];
  double nu = Mm->eliso[im].nu;
  double g0;
  // strain tensor due to positive principal strains
  fillm(0.0, et);
  for (i = 0; i < 3; i++)
  {
    if (peps(i) > 0.0)
      et(i,i) = peps(i);
  }
  lgmatrixtransf(et, t);
  // strain tensor due to negative principal strains
  fillm(0.0, ec);
  for (i = 0; i < 3; i++)
  {
    if (peps(i) < 0.0)
      ec(i,i) = peps(i);
  }
  lgmatrixtransf(ec, t);
  // tensor of dimensionless damage driving force in tension 
  mxm(et, et, tyt);
  g0 = 1.0/(2.0*(1.0+nu));  
  cmulm(g0, tyt);
  glmatrixtransf(tyt, t); // transformation to the principal directions
  // tensor of dimensionless damage driving force in compression 
  mxm(ec, ec, tyc);
  g0 = 1.0/(2.0*(1.0+nu));  
  cmulm(g0, tyc);
  glmatrixtransf(tyc, t);  // transformation to the principal directions
  // transformation to the principal values vector
  for (i = 0; i < 3; i++)
  {   
    pyt(i) = tyt(i,i);
    pyc(i) = tyc(i,i);
  }  
}
/**
  Function computes value of load function for the volumetric damage.

  @param  ipp  - integration point number
  @param  peps - %vector of principal strains
  @param  d    - actual value of damage parameter
  @param  aa   - actual value of material parameter a

  @retval The function returns value of load function for the volumetric damage.
  
  TKo 9.2006
*/
double anisodam::loadfuncvol(long ipp, vector &peps, double d, double aa)
{
  double y, f;
  double tmp1, tmp2;

  y = damdrvforce_vol(ipp, peps);
  if (y <= y0) 
   f = -1.0;      
  else
  {
    // compute (y-y0)^b
    tmp1 = y - y0;
    tmp2 = log(tmp1);
    tmp1 = exp(b*tmp2);
    // compute value of load function
    f = (1.0 - d)*(1.0+aa*tmp1) - 1.0;
  }
  return f;
}


/**
  Function computes value of load function for deviatoric damage

  @param  ipp  - integration point number
  @param  peps - %vector of principal strains
  @param  t    - %matrix of principal directions stored in colomns for given
                 principal strains
  @param  damt - %vector actual values of principal damage parameters for tension
  @param  damc - %vector actual values of principal damage parameters for compression
  @param  aat  - actual values of material parameter a for tension 
  @param  aac  - actual values of material parameter a for compression 
  @param  lft  - %vector of resulting values of load function for tension for each principal direction
  @param  lfc  - %vector of resulting values of load function for compression for each principal direction

  @retval The function returns load function values separately for tension and compression
          via parameters lft and lfc for each principal direction.

  TKo 9.2006
*/
void anisodam::loadfuncdev(long ipp, vector &peps, matrix &t, vector &damt, vector &damc, 
                           double aat, double aac, vector &lft, vector &lfc)
{
  long i;
  double tmp1, tmp2;
  vector pyt(3), pyc(3);

  damdrvforce_dev(ipp, peps, t, pyt, pyc);
  for (i=0; i<3; i++)
  {
    // for tension
    if (pyt(i) <= y0t) 
      lft(i) = -1.0;      
    else
    {
      // compute (yt-y0t)^bt
      tmp1 = pyt(i) - y0t;
      tmp2 = log(tmp1);
      tmp1 = exp(bt*tmp2);
      // compute value of load function for tension in given principal direction i
      lft(i) = (1.0 - damt(i))*(1.0+aat*tmp1) - 1.0;
    }
    // for compression
    if (pyc(i) <= y0c) 
      lfc(i) = -1.0;      
    else
    {
      // compute (yc-y0c)^bc
      tmp1 = pyc(i) - y0c;
      tmp2 = log(tmp1);
      tmp1 = exp(bc*tmp2);
      // compute value of load function for compression in given principal direction i
      lfc(i) = (1.0 - damc(i))*(1.0+aac*tmp1) - 1.0;
    }
  }  
  return;
}


/**
  Function computes increments of damage parameter d. The damage parameter
  is expressed as derivative of damage function with respect to damage parameter d.

  @param ipp - integration point number
  @param y   - driving forces for volumetric damage
  @param dy  - increment of driving forces for volumetric damage
  @param aa  - actual value of material parameter a 
  @param lf  - actual value of load function

  @retval The function returns value of the damage parameter increment.
  
  TKo 9.2006
*/
double anisodam::daminc_vol(long ipp, double y, double dy, double aa, double lf)
{
  double e = Mm->give_actual_ym(ipp);
  double tmp, ddam = 0.0;

  if ((lf > 0.0) && (dy > 0.0))
  {  
    tmp = log(y-y0);
    ddam = aa*b*exp(tmp*(b-1.0));
    tmp = 1.0+a*exp(tmp*b);
    ddam /= e*tmp*tmp;
    ddam *= dy;
  }
  return ddam;
}


/**
  Function computes increments of damage parameters Dt and Dc. The damage parameters
  are expressed as derivatives of damage function with respect to damage parameters .

  @param ipp - integration point number
  @param pyc - %vector of principal driving forces for tension
  @param pyt - %vector of principal driving forces for tension
  @param dyc - %vector of principal driving forces increments for tension
  @param dyt - %vector of principal driving forces increments for compression
  @param aat - actual value of material parameter at 
  @param aac - actual value of material parameter ac 
  @param lft - %vector of actual load function values for tension
  @param lfc - %vector of actual load function values for compression
  @param pdamc - %vector of principal damage parameters increments for tension
  @param pdamt - %vector of principal damage parameters increments for compression
  
  TKo 9.2006
*/
void anisodam::pdaminc_dev(long ipp, vector &pyc, vector &pyt, vector &dyt, vector &dyc, 
                           double aat, double aac, vector &lft, vector &lfc, vector &dpdamt, vector &dpdamc)
{
  long i;
  double e = Mm->give_actual_ym(ipp);
  double tmp, ltmp;

  fillv(0.0, dpdamt);
  fillv(0.0, dpdamc);

  for(i=0; i<3; i++)
  {    
    // tension 
    if ((dyt[i] > 0.0) && (lft[i] > 0.0))
    {
      ltmp = log(pyt(i) - y0t);
      dpdamt(i) = aat*bt*exp(ltmp*(bt-1.0));
      tmp = 1.0+aat*exp(ltmp*bt);
      dpdamt(i) /= e*(1.0+tmp*tmp);
      dpdamt(i) *= dyt(i);
    }
    // compression 
    if ((dyc[i] > 0.0) && (lfc[i] > 0.0))
    {
      ltmp = log(pyc(i) - y0c);
      dpdamc(i) = aac*bc*exp(ltmp*(bc-1.0));
      tmp = 1.0+aac*exp(ltmp*bc);
      dpdamc(i) /= e*(1.0+tmp*tmp);
      dpdamc(i) *= dyc(i);
    }
  }
}



/**
  Function computes value of damage parameter d. The damage parameter
  is derived from the damage function.

  @param ipp - integration point number
  @param y   - driving forces for volumetric damage
  @param dy  - increment of driving forces for volumetric damage
  @param aa  - actual value of material parameter a 
  @param lf  - actual value of load function

  @retval The function returns value of the volumetric damage parameter.
  
  TKo 9.2006
*/
double anisodam::dam_vol(long ipp, double y, double dy, double aa, double lf)
{
  double tmp, dam = 0.0;

  if ((lf > 0.0) && (dy > 0.0))
  {  
    tmp = log(y-y0);
    dam = aa*exp(tmp*b);
    tmp = 1.0/(1.0+dam);
    dam = 1.0-tmp;
  }
  return dam;
}


/**
  Function computes values of principal damage parameters Dt and Dc. The damage parameters
  are derived from the damage function.

  @param ipp - integration point number
  @param pyc - %vector of principal driving forces for tension
  @param pyt - %vector of principal driving forces for tension
  @param dyc - %vector of principal driving forces increments for tension
  @param dyt - %vector of principal driving forces increments for compression
  @param aat - actual value of material parameter at 
  @param aac - actual value of material parameter ac 
  @param lft - %vector of actual load function values for tension