creep_b3.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_b3.h"

b3mat::b3mat (void)
{
  q1=q2=q3=q4=q5=0.0;

  e28=30.0e9; //Pa
  fc=35.8e6; //Pa
  ft=1.5e6; //Pa
  ft_ratio = 1.0e-5;
  alpha = 12.0e-6; //K^-1
  wc=0.43;
  sc=3.4;
  gc=1.98;
  cs=305.0;  //kg*m^-3
  a1=1.05;
  a2=1.2;
  kd=0.15;   //m
  type_e=0;
  type_h=0;
  type_temp=0;

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

  timeMax=(Mp->timecon.endtime ())/86400.0;//in days
      
  flag_drshr = flag_shr = flag_temp = 0;

  retTime=NULL;

  eps_ainf = 0.0;
}

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

/**
   function reads material parameters

   @param in - input file

   TKr, 3.6.2005
*/

void b3mat::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 %lf %lf %lf %lf %lf %lf %ld %ld", &fc, &ft, &ft_ratio, &alpha, &wc, &sc, &gc, &cs, &a1, &a2, &ks, &kd, &type_h, &type_temp);
  
  //input in Pa
  fc = fc/6.89476*1.0e-3;//to psi
  //input in kg*m^-3
  cs = cs/16.03;//to lb*ft^-3
  //input in m
  kd=kd*1000.0/25.4;//to inch
  
  //environmental humidity and temperature
  if (type_h == 0){xfscanf (in,"%lf ",&hum_env); }
  if (type_temp == 0){xfscanf (in,"%lf ",&temp_env); }

 xfscanf (in,"%lf %lf %ld %ld %ld %d %d %d",
	 &tb_time,&th_time,&napproxtime,&nRetTime,&type_rt,&flag_drshr,&flag_shr,&flag_temp);

 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]);
 }

 xfscanf (in,"%lf",&eps_ainf);

}




/**
   function reads material parameters

   @param out - output file

   TKr, 9.5.2007
*/

void b3mat::print (FILE *out)
{
  long i;

  fprintf(out," %ld",type_e);
  
  if (type_e == 1){fprintf (out," %e ",e28);}
  
  fprintf (out," %e %e %e %e %e %e %e %e %e %e %e %e %ld %ld",fc,ft,ft_ratio,alpha,wc,sc,gc,cs,a1,a2,ks,kd,type_h,type_temp);
  
  if (type_h == 0){fprintf (out," %lf ",hum_env); }
  if (type_temp == 0){fprintf (out," %lf ",temp_env); }

  fprintf (out," %e %e %ld %ld %ld %d %d %d",
	 tb_time,th_time,napproxtime,nRetTime,type_rt,flag_drshr,flag_shr,flag_temp);

 if (type_rt==1){//retardation times (ages) [days]
   for (i=0;i<nRetTime;i++)
     fprintf (out," %e",retTime[i]);
 }
 fprintf (out," %e",eps_ainf);
}

/**
   function converts times to ages

   @param ipp - number of integration point

   TKr, 3.6.2005
*/
void b3mat::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 17.7.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;

  if (tb_age <= 0.0){
    fprintf (stderr,"\n Age of concrete must be greater than zero!!!");
    fprintf (stderr,"\n Age of concrete = %e, Actual time = %e, Time of end of concrete casting = %e; (Age of concrete = Actual time - Time of end of concrete casting)",tb_age,actualtime,tb_time);
    fprintf (stderr,"\n In function give_ages (file %s, line %d).\n\n\n\n\n\n\n\n",__FILE__,__LINE__);
    exit(0);
  }
  if (th_age <= 0.0){
    fprintf (stderr,"\n Age of concrete when drying begins must be greater than zero!!!");
    fprintf (stderr,"\n Age of concrete when drying begins = %e, Actual time = %e, Time when drying begins = %e; (Age of concrete when drying begins = Actual time - Time when drying begins).",th_age,actualtime,th_time);
    fprintf (stderr,"\n In function give_ages (file %s, line %d).\n\n\n\n\n\n\n\n",__FILE__,__LINE__);
    exit(1);
  }

  if (th_time <= tb_time){
    fprintf (stderr,"\n Time when drying begins must be greater than Time of end of concrete casting!!!");
    fprintf (stderr,"\n Time when drying begins = %e   Time of end of concrete casting = %e.",th_time,tb_time);
    fprintf (stderr,"\n In function give_ages (file %s, line %d).\n\n\n\n\n\n\n\n",__FILE__,__LINE__);
    exit(2);
  }

  t_0 = th_time - tb_time;
}


/**
   function returns number of components of other array

   @param ipp - number of integration point

   TKr, 3.6.2005
*/

long b3mat::give_nceqother (long ipp)
{
  long nc,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;           //stress-induced shrinkage-irreversible strains increments
  nceqother = nceqother + nc;           //total irreversible strains (creep-irreversible strains + shrinkage-irreversible strains + stress-induced shrinkage-irreversible strains)
  nceqother = nceqother + 1;            //previous free shrinkage
  
  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 b3mat::initvalues (long ipp, long ido)
{
  long nc;

  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 b3mat::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 b3mat::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 b3mat::store_emu(vector e_mu,long n_ret_times,long ipp,long ido)
{
  long i;
  long nc;

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

  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 b3mat::give_emu(vector &e_mu,long n_ret_times,long ipp,long ido)
{
  long i;
  long nc;

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

  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 b3mat::give_hidden_strains(matrix &gamma_mu,long ipp,long ido)
{
  long i,j,ii;
  long n_ret_times;
  long nc;

  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 b3mat::store_hidden_strains(matrix gamma_mu,long ipp,long ido)
{
  long i,j,ii;
  long n_ret_times;
  long nc;

  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 other

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

   TKr, 3.6.2005
*/
void b3mat::give_stresses(vector &sigma,long ipp,long ido)
{
  long i;
  long nc;

  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 other

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

   TKr, 3.6.2005
*/
void b3mat::store_stresses(vector sigma,long ipp,long ido)
{
  long i;
  long nc;

  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 other

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

   TKr, 3.6.2005
*/
void b3mat::give_dstresses(vector &dsigma,long ipp,long ido)
{
  long i;
  long nc;

  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 other

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

   TKr, 3.6.2005
*/
void b3mat::store_dstresses(vector dsigma,long ipp,long ido)
{
  long i;
  long nc;

  nc=Mm->ip[ipp].ncompstr;
  
  for (i=0;i<nc;i++){
    Mm->ip[ipp].eqother[ido+(nc+nc+nc)+i] = dsigma[i];
  }
}