creep_rspec.cpp

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

/** 
    Continuos retardation spectrum for solidification theory of concrete creep for B3 model
    according to Bazant and Prasannan and V. Smilauer
*/

#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_rspec.h"



rspecmat::rspecmat (void)
{
  n_param = 0.0;
  m_param = 0.0;

  q1=q2=q3=q4=q5=0.0;
  C_const = 0.0;

  min_ret_time = 1.e-5;//minimum retardation time used in calculation [days], truncated values will be hidden in C_{const}, should be about 1.0e-16 for exact q1 determination, otherwise can be like 1.e-3

  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 = 10;
  type_rt = 0;
  e0 = 0.0;
  actualtime = 0.0;
  dt = 0.0;

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

  eps_ainf = 0.0;
}

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

/**
   function reads material parameters

   @param in - input file

   TKr, 3.6.2005
*/

void rspecmat::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/0.0254;//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 rspecmat::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 rspecmat::give_ages (double &tb_age_dt,double &tb_age,double &tl_age,double &th_age,double &dt,double &maxtime,long &napptime,long ipp)
{
  //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__);
    //print (Out);//debug??!!
    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 rspecmat::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 rspecmat::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 rspecmat::give_nret_time (void)
{
  return(nRetTime);
}

double rspecmat::give_L(double q3, double q, double tau){

  double L;

  L = -2.*n_param*n_param*pow(3.*tau, 2.*n_param-3.)*(n_param-1.-pow(3.*tau,n_param)) / pow((1.+pow(3.*tau,n_param)),3.);
  L += (n_param*(n_param-2.)*pow(3*tau, n_param-3.)*(n_param-1.-pow(3.*tau,n_param))-n_param*n_param*pow(3.*tau,2.*n_param-3.)) / ((1.+pow(3.*tau,n_param))*(1.+pow(3.*tau,n_param)));
  L *= q3*pow(3.*tau,3.)/2.;
  L *= log (10.) * q;//now L is A_mu = 1/E_mu
  
  return L;
}



/**
   function computes creep compliacne function J(t,t') [1/Pa] in a material point B3 Bazant's model
   dependent on temerature and humidity changes

   @param t0        - age when drying begins [days]
   @param tl        - age at loading [days]
   @param t         - time representing age of concrete [days]
   @param ipp       - number of integration point
  
   TKr, 4.9.2006
*/
//jeste doplnit efekt pro nelinearni creep??!!
double rspecmat::give_J_E_mu(vector &e_mu,double t0, double tl, double t, long ipp,long ido){
  
  double jt,L,ret_time_hlp;
  double n,cd,ac,ag,m,z,r,qf,q,eps_shinf,tau,ht,htl,st,stl;
  double hum,temp,hum_prev,temp_prev;
  double kt;
  long nc;

  //added 22.6.2007 correction:
  t0 = t_0;

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

  if ( type_h == 0){
    hum = 0.0;
    hum_prev = 0.0;
  }
  else{
    // actual humidity
    hum=Mm->moist[ipp];
    // humidity from previous time step
    hum_prev = Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc)];
  }
  if ( type_temp == 0){
    temp = 0.0;
    temp_prev = 0.0;
  }
  else{
    // actual temperature
    temp=Mm->tempr[ipp];
    // temperature from previous time step
    temp_prev = Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc)+1];
  }

  if (t < t0){
    fprintf (stderr,"\n Age of concrete is lower than age when dyring begins!!!");
    fprintf (stderr,"\n Age of concrete = %e,   Age when dyring begins = %e",t,t0);
    fprintf (stderr,"\n In function give_J_E_mu (file %s, line %d).\n",__FILE__,__LINE__);
    abort();
  }
  if (t < tl){
    fprintf (stderr,"\n Age of concrete is lower than age at loading!!!");
    fprintf (stderr,"\n Age of concrete = %e,   Age at loading = %e",t,tl);
    fprintf (stderr,"\n In function give_J_E_mu (file %s, line %d).\n",__FILE__,__LINE__);
    abort();
  }

  //basic creep
  ac=sc+gc;   ag=ac/gc;   m_param = 0.5; m = m_param; //m=0.28+1.0/fc/fc;
  n_param = 0.1;
  n = n_param;

  if (type_e == 1)
    q1=600000.0/e28;//measured
  else
    q1=600000.0/57000.0/sqrt(fc);//empirical
  r=1.7*pow(tl,0.12)+8.0;
  qf=1.0/(0.086*pow(tl,(2.0/9.0))+1.21*pow(tl,(4.0/9.0)));
  z = 1.0/pow(tl,m)*log(1.0 + pow((t - tl),n));
  q = qf/pow((1.0 + pow((qf/z),r)),(1.0/r));
  q2=451.1*sqrt(cs)/pow(fc,0.9);
  q3=0.29*pow(wc,4.0)*q2;
  q4=0.14/pow(ac,0.7);
  
  //drying creep
  cd = 0.0;
  //humidity effect
  if (type_h == 0){
    //constant huminidy
    //kt = 190.8/pow(t0,0.08)*fc;//this is wrong this is from literature (t0 = age when drying begins)
    kt = 190.8/pow(t,0.08)/pow(fc,0.25);//this is according to creepb.cpp (Fajman)
    tau = kt*ks*ks*kd*kd;
    st = tanh(sqrt((t - t0)/tau));
    stl = tanh(sqrt((tl - t0)/tau));
    ht = 1.0 - (1.0 - hum_env)*st;
    htl = 1.0 - (1.0 - hum_env)*stl;
    eps_shinf = a1*a2*(26.0*pow((wc*cs),2.1)/pow(fc,0.28)+270.);
    q5 = 7.57e5/fc/pow(eps_shinf,0.6);
    cd = q5*sqrt(exp(-8.0*ht) - exp(-8.0*htl)); 

    if (tl < t0)
      cd = 0.0;    
  }
  /*   if (type_h == 1){
       //changing humidity
       //stress-induced effect is included as strain in function: give_deps_stressinduced, this is old:
       dhum = hum-hum_prev; 
       eps_shinf = a1*a2*(26.0*pow((wc*cs),2.1)/pow(fc,0.28)+270.);
       q5 = 0.35/fc;
       cd = fabs(q5*eps_shinf*3.0*hum*hum*dhum);
       }
       
       //temperature effect
       if ( type_temp == 0 ){
       }
       else {
       //changing temperature
       dtemp = temp - temp_prev;   
       q5= 1.5/fc;
       cd = cd+fabs(q5*alpha*dtemp);
       }
  */
  if (flag_drshr == 0)
    cd = 0.0;

  //units:
  q1 = q1/6.89476*1.0e-3*1.0e-6;//units 1e-6*psi^-1 -> Pa^-1
  q2 = q2/6.89476*1.0e-3*1.0e-6;//units 1e-6*psi^-1 -> Pa^-1
  q3 = q3/6.89476*1.0e-3*1.0e-6;//units 1e-6*psi^-1 -> Pa^-1
  q4 = q4/6.89476*1.0e-3*1.0e-6;//units 1e-6*psi^-1 -> Pa^-1
  q5 = q5/6.89476*1.0e-3*1.0e-6;//units 1e-6*psi^-1 -> Pa^-1

  q = log10(retTime[1]/retTime[0]);//log spacing of consequent retardation times  
  jt = 0.;
  C_const = 0;
  
  //calculate degenerated (already crept) Kelvin series from min_ret_time to 10^-20 days
  ret_time_hlp = retTime[0];
  while (ret_time_hlp >= 1.e-20){
    ret_time_hlp /= pow(10.,q);//calculate next retardation time which is smaller in log scale
    C_const += give_L(q3, q, ret_time_hlp);
  }

  for(int i=0; i<nRetTime; i++){
    L = give_L(q3, q, retTime[i]);
    e_mu[i] = 1./L;
    jt+=L*(1.-exp(-(t-tl)/retTime[i]));
  }

  jt += q1 + C_const;

  //returns only e_mu
  return (jt);
}

//q1
double rspecmat::give_q1()
{
  return(q1);
}

//C_const
double rspecmat::give_C_const()
{
  return(C_const);
}


//q4
double rspecmat::give_q4()
{ 
  return(q4);
}


//nonlin_func
double rspecmat::give_nonlin_func()
{
  double nonlin_func;
  
  nonlin_func = 1.0;

  return(nonlin_func);
}


//inv_v
double rspecmat::give_inv_v(double time_mid){
  
  double inv_v;

  inv_v = q2/q3 * pow(1.0/time_mid, m_param) + 1.;

  return(inv_v);
}

/**
   function computes retardation times

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

*/
void rspecmat::give_rettimes (vector &rettimes,long n_ret_times,long ipp)
{
  long i;

  //retardation times on a spectrum of non-aging log-power law from Laplace transformation
  rettimes.a[0] = retTime[0] = min_ret_time;
  
  for (i=1;i<nRetTime;i++){
    rettimes.a[i] = retTime[i] = retTime[0] * pow(10.,i);//less spacing does not make any significant improvement 
  }
}


/**
   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 rspecmat::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 rspecmat::give_emu(vector &e_mu,long n_ret_times,long ipp,long ido)
{