creep_dpl.cpp

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

#include "global.h"
#include "math.h"
#include "globmat.h"
#include "genfile.h"
#include "adaptivity.h"
#include "intpoints.h"
#include "stdlib.h"
#include "elastisomat.h"
#include "creep_dpl.h"

dplmat::dplmat (void)
{
  e28=30.0e9; //Pa
  fc=35.8e6; //Pa
  ro = 2500.0; //kg*m-3
  wc=0.43;
  sc=3.4;
  gc=1.98;
  a1 = 1.05;

  tb_time = 0.0;
  th_time = 0.0;
  napproxtime = 0;
  nRetTime = 0;
  type_rt = 0;
  e0 = 0.0;
  actualtime = 0.0;
  dt = 0.0;
  nc = 0;

  timeMax=(Mp->timecon.endtime ())/86400.0;//in days

  retTime=NULL;
}

dplmat::~dplmat (void)
{
  delete [] retTime;
}

/**
   function reads material parameters

   @param in - input file

   TKr, 3.6.2005
*/

void dplmat::read (XFILE *in)
{
  long i;

  xfscanf(in,"%ld",&type_e);
  
  if (type_e == 1){xfscanf (in,"%lf ",&e28);}
  //input in Pa
  e0 = 1.5*e28/6.89476*1.0e-3;//to psi
  e28 = e28/6.89476*1.0e-3;//to psi
  
  xfscanf (in,"%lf %lf %lf %lf %lf %lf",&fc,&ro,&wc,&sc,&gc,&a1);
  
  //input in Pa
  fc = fc/6.89476*1.0e-3;//to psi
  //input in kg*m^-3
  ro = ro/16.03;//to lb*ft^-3
  
  xfscanf (in,"%lf %lf %ld %ld %ld",
	  &tb_time,&th_time,&napproxtime,&nRetTime,&type_rt);
  
  tb_time = tb_time/86400.0;//days
  th_time = th_time/86400.0;//days
  
  retTime = new double [nRetTime];//vector of ret. times
  for (i=0;i<nRetTime;i++){
    retTime[i]=0.0;
  }
  
  if (type_rt==1){//reading of retardation times (ages) [days]
    for (i=0;i<nRetTime;i++)
      xfscanf (in,"%lf",&retTime[i]);
  }
}


/**
   function converts times to ages

   @param ipp - number of integration point

   TKr, 3.6.2005
*/
void dplmat::give_ages (double &tb_age_dt,double &tb_age,double &tl_age,double &th_age,double &dt,double &maxtime,long &napptime,long ipp)
{
  //actuall times in days
  /*   //old  
       dt = Mp->timecon.actualforwtimeincr ()/86400.0;//days
       actualtime = Mp->timecon.actualtime ()/86400.0;//days
       th_age = th_time;
       tl_age = actualtime + dt/2.0;
       tb_age = actualtime;
       napptime = give_napproxtime(ipp);
       maxtime = Mp->timecon.endtime ()/86400.0;//days
       tb_age_dt = tb_age + dt;
  */
  
  //new 22.8.2006
  dt = Mp->timecon.actualforwtimeincr ()/86400.0;//days
  actualtime = Mp->timecon.actualtime ()/86400.0;//days
  th_age = actualtime - th_time;
  tl_age = actualtime - tb_time + dt/2.0;
  tb_age = actualtime - tb_time;
  napptime = give_napproxtime(ipp);
  //oprava 6.11.2006
  //maxtime = Mp->timecon.endtime ()/86400.0 - Mp->timecon.starttime ()/86400.0;//days
  maxtime = Mp->timecon.endtime ()/86400.0;//days
  tb_age_dt = tb_age + dt;
}



/**
   function returns number of components of other array

   @param ipp - number of integration point

   TKr, 3.6.2005
*/

long dplmat::give_nceqother (long ipp)
{
  long nceqother;

  nc=Mm->ip[ipp].ncompstr;

  nceqother = nc + nc + nc + nc;        //total strains, total stresses, total strain increments, total stress increments
  nceqother = nceqother + 2 + 2;        //total (actual) humidity, total (actual) temperature, humidity increment, temperature increment
  //  nceqother = nceqother + nRetTime;     //retardation times
  nceqother = nceqother + nRetTime;     //coefficients of Dirichlet series
  nceqother = nceqother + nc*nRetTime;  //hidden strains increments
  nceqother = nceqother + nc + nc;      //creep-irreversible strains increments, shrinkage-irreversible strains increments
  nceqother = nceqother + nc;           //irreversible strains (creep-irreversible strains + shrinkage-irreversible strains)

  return(nceqother);
}



/**
   Function initializes eqother array with initial temperature.
   Actual value of array Mm->tempr[ipp] is taken as initial temperature.
   
   @param ipp - integration point pointer
   @param ido - index of internal variables for given material in the ipp other array
   
   21.6.2005, TKo+TKr
*/
void dplmat::initvalues (long ipp, long ido)
{
  nc=Mm->ip[ipp].ncompstr;
  if(Mm->moist != NULL)
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc)]=Mm->moist[ipp];
  if(Mm->tempr != NULL)
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+1)]=Mm->tempr[ipp];
}




/**
   function returns number of retardation times

   TKr, 3.6.2005
*/

long dplmat::give_nret_time (void)
{
  return(nRetTime);
}




/**
   function computes retardation times

   @param rettimes - vector of ret. times
   @param ipp      - number of integration point

   TKr, 3.6.2005
*/
void dplmat::give_rettimes (vector &rettimes,long n_ret_times,long ipp)
{
  long i;
  double m,mm;

  if (n_ret_times != nRetTime){//check of number of ret. times
    fprintf (stderr,"\n Wrong number of ret. times in function give_rettimes (file %s, line %d).\n",__FILE__,__LINE__);
    abort();
  }

  if (type_rt==1){
  }
  else
    {
      m = log10(2.0*timeMax);
      
      retTime[0]=1.0e-9;
      retTime[1]=1.0;
      mm=1./(nRetTime-2);
      
      for (i=2;i<nRetTime-1;i++){
	retTime[i]=pow(10.,(i-1)*mm*m);
      }
      retTime[nRetTime -1]=2.0*timeMax;
    }
  
  for (i=0;i<nRetTime;i++){
    rettimes.a[i]=retTime[i];
  }
}



/**
   function stores E_mu stifnesses of Kelvin chains

   @param e_mu      - vector of stifnesses
   @param nrettime - number of retardation times
   @param ipp      - number of integration point

   TKr, 3.6.2005
*/
void dplmat::store_emu(vector e_mu,long n_ret_times,long ipp,long ido)
{
  long i;

  if (n_ret_times != nRetTime){//check of number of ret. times
    fprintf (stderr,"\n Wrong number of ret. times in function store_emu (file %s, line %d).\n",__FILE__,__LINE__);
    abort();
  }
  
  for (i=0;i<nRetTime;i++){
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2)+ i]  = e_mu.a[i];
  }
}

/**
   function returns E_mu stifnesses of Kelvin chains

   @param e_mu      - vector of stifnesses
   @param nrettime - number of retardation times
   @param ipp      - number of integration point

   TKr, 3.6.2005
*/
void dplmat::give_emu(vector &e_mu,long n_ret_times,long ipp,long ido)
{
  long i;

  if (n_ret_times != nRetTime){//check of number of ret. times
    fprintf (stderr,"\n Wrong number of ret. times in function store_emu (file %s, line %d).\n",__FILE__,__LINE__);
    abort();
  }
  
  for (i=0;i<nRetTime;i++){
    e_mu.a[i] = Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2)+ i];
  }
}



/**
   function returns Gamma_mu hidden strains of Kelvin chains

   @param gamma_mu - vector of stifnesses
   @param nrettime - number of retardation times
   @param nc       - number of strain components
   @param ipp      - number of integration point

   TKr, 3.6.2005
*/
void dplmat::give_hidden_strains(matrix &gamma_mu,long ipp,long ido)
{
  long i,j,ii;
  long n_ret_times;
  
  nc=Mm->ip[ipp].ncompstr;
  n_ret_times = give_nret_time ();

  ii = 0;
  for (i=0;i<n_ret_times;i++){
    for (j=0;j<nc;j++){
      gamma_mu[j][i]=Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times)+ ii];
    ii++;
    }
  }
}



/**
   function stores Gamma_mu hidden strains of Kelvin chains

   @param gamma_mu - vector of stifnesses
   @param nrettime - number of retardation times
   @param nc       - number of strain components
   @param ipp      - number of integration point

   TKr, 3.6.2005
*/
void dplmat::store_hidden_strains(matrix gamma_mu,long ipp,long ido)
{
  long i,j,ii;
  long n_ret_times;
  
  nc=Mm->ip[ipp].ncompstr;
  n_ret_times = give_nret_time ();

  ii = 0;
  for (i=0;i<n_ret_times;i++){
    for (j=0;j<nc;j++){
      Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times)+ ii] = gamma_mu[j][i];
    ii++;
    }
  }
}





/**
   function returns total stresses from ip

   @param gamma_mu - vector of total stresses
   @param nc       - number of components
   @param ipp      - number of integration point

   TKr, 3.6.2005
*/
void dplmat::give_stresses(vector &sigma,long ipp,long ido)
{
  long i;
  
  nc=Mm->ip[ipp].ncompstr;
  
  for (i=0;i<nc;i++){
    sigma[i] = Mm->ip[ipp].eqother[ido+(nc)+i];
  }
}



/**
   function stores total stresses into ip

   @param gamma_mu - vector of total stresses
   @param nc       - number of components
   @param ipp      - number of integration point

   TKr, 3.6.2005
*/
void dplmat::store_stresses(vector sigma,long ipp,long ido)
{
  long i;
  
  nc=Mm->ip[ipp].ncompstr;
  
  for (i=0;i<nc;i++){
    Mm->ip[ipp].eqother[ido+(nc)+i] = sigma[i];
  }
}


/**
   function returns increments of total stresses from ip

   @param dsigma - vector of increments of total stresses
   @param nc     - number of components
   @param ipp    - number of integration point

   TKr, 3.6.2005
*/
void dplmat::give_dstresses(vector &dsigma,long ipp,long ido)
{
  long i;
  
  nc=Mm->ip[ipp].ncompstr;
  
  for (i=0;i<nc;i++){
    dsigma[i] = Mm->ip[ipp].eqother[ido+(nc+nc+nc)+i];
  }
}



/**
   function stores increments of total stresses into ip

   @param sigma - vector of increments of total stresses
   @param nc    - number of components
   @param ipp   - number of integration point

   TKr, 3.6.2005
*/
void dplmat::store_dstresses(vector dsigma,long ipp,long ido)
{
  long i;
  
  nc=Mm->ip[ipp].ncompstr;
  
  for (i=0;i<nc;i++){
    Mm->ip[ipp].eqother[ido+(nc+nc+nc)+i] = dsigma[i];
  }
}




/**
   function stores increments of total stresses into ip and adds them to total stresses

   @param sigma - vector of increments of total stresses
   @param nc    - number of components
   @param ipp   - number of integration point

   TKr, 3.6.2005
*/
void dplmat::store_dstresses_add(vector dsigma,long ipp,long ido)
{
  long i;
  
  nc=Mm->ip[ipp].ncompstr;
  
  for (i=0;i<nc;i++){
    Mm->ip[ipp].eqother[ido+(nc+nc+nc)+i] = dsigma[i];
    Mm->ip[ipp].eqother[ido+(nc)+i] = Mm->ip[ipp].eqother[ido+(nc)+i] + dsigma[i];
    Mm->ip[ipp].stress[i] = Mm->ip[ipp].eqother[ido+(nc)+i];
  }
}




/**
   function returns total strains from ip

   @param eps - vector of total strains
   @param nc  - number of components
   @param ipp - number of integration point