creep_b3.cpp

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

    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);
  
  c0=q2*q+q3*log(1.0+pow((t-tl),n))+q4*log(t/tl);

  //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;
  
  //jt=(q1+c0+cd)/6.89476*1.0e-3*1.0e-6;//units 1e-6*psi^-1 -> Pa^-1//old
  //new - only time function (without q1)
  q1 = q1/6.89476*1.0e-3*1.0e-6;//units 1e-6*psi^-1 -> Pa^-1
  jt=(c0+cd)/6.89476*1.0e-3*1.0e-6;//units 1e-6*psi^-1 -> Pa^-1
  
  return(jt);
}


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


/**
   function computes time increment of irreversible free strains in a material point by B3 Bazant's model
   from temerature and humidity changes (drying shrikage, shrinkage, stress induced strain etc.)
   free thermal strain is not included!!      

   @param t0        - age when drying begins [days]
   @param t_dt      - time representing age of concrete + dt [days]
   @param t         - time representing age of concrete [days]
   @param ipp       - number of integration point

   TKr, 3.6.2005
*/
void b3mat::give_deps_free (vector &deps_sh, double t0, double t_dt, double t, double dt, long ipp,long im,long ido)
{
  double eps_shinf,kh,kt,tau,st,dhum,dtemp;
  double hum,temp,hum_prev,temp_prev;
  double deps,deps_h,deps_a,deps_h_old,deps_t,e_607,e_t0_tau,e_t0,e_actual;
  long nc;

  nc=Mm->ip[ipp].ncompstr;  
  
  //added 22.6.2007 correction:
  t0 = t_0;

  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 = %10.15e,   Age when dyring begins = %10.15e",t,t0);
    fprintf (stderr,"\n In function give_deps_free (file %s, line %d).\n",__FILE__,__LINE__);
    abort();
  }

  //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));
    //derivative of st with respect to time
    //dst = tanh(sqrt((t_dt - t0)/tau)) - tanh(sqrt((t - t0)/tau));
    eps_shinf = a1*a2*(26.0*pow((wc*cs),2.1)/pow(fc,0.28)+270.);
    //time dependence of ultimate shrinkage = effect of aging
    e_607 = e28*sqrt(607.0/(4+0.85*607.0));// empirical formulation
    e_t0_tau = e28*sqrt((t0+tau)/(4+0.85*(t0+tau)));// empirical formulation
    eps_shinf = eps_shinf*e_607/e_t0_tau;
    

    if(hum_env<=0.98)
      kh=(1.0-pow(hum_env,3.0));
    else 
      kh= (-0.2 - (1.0-pow(0.98,3.0)))/0.02*(hum_env - 0.98);
    if(hum_env>=1.0)
      kh = -0.2;

    deps_h = -eps_shinf*kh*st; //shrinkage strain

    deps_h_old = give_shrinkage(ipp,ido);
    store_shrinkage(deps_h,ipp,ido);
    deps_h = deps_h - deps_h_old; //shrinkage strain increment

    deps_h = deps_h*1.e-6;//
  }    
  else{
    //changing humidity
    dhum = hum - hum_prev;
    eps_shinf = a1*a2*(26.0*pow((wc*cs),2.1)/pow(fc,0.28)+270.);
    eps_shinf = eps_shinf*1.e-6;//units
    //eps_shinf = 0.0008;//new 24.2.2006 corresponds to Pa,m,N

    //time dependence of ultimate shrinkage
    e_t0 = e28*sqrt((t0)/(4+0.85*(t0)));// empirical formulation
    e_actual = give_actual_ym (ipp,im,ido)/6.89476*1.0e-3;//to psi
    //eps_shinf = eps_shinf*e_t0/e_actual;//zatim??!!

    deps_h = eps_shinf*3.0*hum*hum*dhum; //shrinkage strain increment

    //simple empirical formula for autogenous shrinkage
    //added total autogenous shrinkage
    deps_a = -eps_ainf*((1.0-exp(-.125*t_dt*24.0))-(1.0-exp(-.125*t*24.0)));//time in hours
    
    deps_h = deps_h + deps_a;
  }
  

  //temperature effect
  dtemp = temp - temp_prev;  
  deps_t = 0.0;//deps_t = alpha*dtemp;//free thermal strain is not included now
  
  if (flag_shr == 0)
    deps_h = 0.0;
  
  if (flag_temp == 0)
    deps_t = 0.0;
  
  deps = deps_h + deps_t;
   
  fillv(0.0,deps_sh);
  ss=Mm->ip[ipp].ssst;
  
  //only uniaxial
  if(ss==bar){
    deps_sh[0]=deps;
  }
  else{
    deps_sh[0]=deps;
    deps_sh[1]=deps;
  }
  if(ss==planestrain) deps_sh[3]=deps;
  if(ss==planestress) deps_sh[3]=deps;
  if(ss==axisymm) deps_sh[2]=deps;
  if(ss==spacestress) deps_sh[2]=deps;

}



/**
   function computes time increment of irreversible stress induced strains in a material point B3 Bazant's model
   from temerature and humidity changes (drying shrikage, shrinkage, stress induced strain etc.)
   
   @param t0        - age when drying begins [days]
   @param sigma     - actual uniaxial stress [Pa]
   @param t_dt      - time representing age of concrete + dt [days]
   @param t         - time representing age of concrete [days]
   @param ipp       - number of integration point

   TKr, 3.6.2005
*/
void b3mat::give_deps_stressinduced (vector &deps_ss, double t0, double t_dt, double t, vector sigma, long ipp,long im,long ido)
{
  double eps_shinf,deps,deps_cs,deps_ts,r,gh,dhum,rho,gt,dtemp,k;
  double hum,temp,hum_prev,temp_prev,e_t0,e_actual;
  long i,nc;

  nc=Mm->ip[ipp].ncompstr;
  
  //added 22.6.2007 correction:
  t0 = t_0;
  
  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 In function deps_stressinduced (file %s, line %d).\n",__FILE__,__LINE__);
    abort();
  }

  eps_shinf = a1*a2*(26.0*pow((wc*cs),2.1)/pow(fc,0.28)+270.);
  eps_shinf = eps_shinf*1.e-6;//units
  //eps_shinf = 0.0008;//new 24.2.2006 corresponds to Pa,m,N

  //time dependence of ultimate shrinkage
  e_t0 = e28*sqrt((t0)/(4+0.85*(t0)));// empirical formulation
  e_actual = give_actual_ym (ipp,im,ido)/6.89476*1.0e-3;//to psi
  //eps_shinf = eps_shinf*e_t0/e_actual;//zatim??!!

  k = -eps_shinf;
  r = 0.3*ft*1.e-6;//in MPa^-1
  rho = 1.5*ft*1.e-6;//in MPa^-1

  dhum = hum-hum_prev;
  dtemp = temp - temp_prev;

  if(dhum >= 0.0)
    gh = 1.0;
  else
    gh = -1.0;

  if(dtemp >= 0.0)
    gt = 1.0;
  else
    gt = -1.0;

  
  fillv(0.0,deps_ss);
  ss=Mm->ip[ipp].ssst;  

  for(i=0;i<nc;i++){
    if (type_h == 0)
      deps_cs = 0.0;
    else
      deps_cs = -k*r*gh*sigma[i]*dhum;
    
    if (type_temp == 0)
      deps_ts = 0.0;
    else
      deps_ts = -alpha*rho*gt*sigma[i]*dtemp;
    
    deps = (deps_cs + deps_ts)*1.e-12;//units
      
    if (flag_drshr == 0)
      deps = 0.0;
    
    deps_ss[i]=deps;
  }

  //only uniaxial
  if(ss==planestress){
    deps_ss[2]=0.0;
  }
  if(ss==planestrain){
    deps_ss[2]=0.0;
  }
  if(ss==axisymm){
    deps_ss[3]=0.0;
  }
  if(ss==spacestress){
    deps_ss[3]=0.0;
    deps_ss[4]=0.0;
    deps_ss[5]=0.0;
  }

}


/**
  Function stores increments of irreversible strains into total irreversible strains.

  @param ipp   - integration point number in the mechmat ip array.
  @param ido   - index of the first internal variable for given material in the ipp other array
  @param deps  - %vector of increments of irreversible strains which to be stored into total irr. strains
*/
void b3mat::storeirrstrains (vector deps,long ipp, long ido)
{
  long i;
  long n_ret_times;
  long nc;

  nc=Mm->ip[ipp].ncompstr;
  
  n_ret_times = give_nret_time ();

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



/**
  Function returns total irreversible strains.

  @param ipp   - integration point number in the mechmat ip array.
  @param ido   - index of the first internal variable for given material in the ipp other array
  @param epscr  - %vector of total irreversible strains
*/
void b3mat::giveirrstrains (long ipp, long ido, vector &epscr)
{
  long i;
  long n_ret_times;
  long nc;
  
  nc=Mm->ip[ipp].ncompstr;
  
  n_ret_times = give_nret_time ();
  
  for (i=0;i<epscr.n;i++){
    epscr[i] = Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times+nc*n_ret_times+nc+nc+nc)+i];
  }
}


/**
   function returns other value from other
   @param compother - number of other components
   @param ipp - first integration point on element
*/

double b3mat::get_othervalue(long compother,long ipp)
{
  double other;

  //dodelat??!!
  switch (compother){ 
  case 0:{//eps_x
    other = 0.0;
      break;
  }
  default:{
    fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
  }
  }
  return (other);

}

/**
     function prints name of other value
     @param out - output file
     @param compother - number of other components
*/
void b3mat::print_othervalue_name(FILE *out,long compother)
{
  //dodelat??!!
  switch (compother){
  case 0:{//eps_x
    fprintf (out,"eps_x ()");
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
  }
  }
}



/** 
   function returns actual Young's modulus
    
   @param ipp - index of integration point
   @param im  - index of
   @param ido - index in array other  

*/ 
double b3mat::give_actual_ym (long ipp,long im,long ido) 
{ 
  long imat;
  double q,e_0,e_t;

  imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
  //creep coefficient
  q = Mm->creep_matstiffchange (ipp,im,ido);
  //asympotic elastic modulus 
  e_0=Mm->eliso[imat].e;
  e_t = q*e_0;
  return (e_t);
} 


/** 
   function returns actual tensile strenght
    
   @param ipp - index of integration point
   @param im  - index of
   @param ido - index in array other  

*/ 
double b3mat::give_actual_ft (long ipp,long im,long ido) 
{ 
  /* double fft,hour;
     
  hour = (Mp->timecon.endtime () - Mp->time)/3600.0;
  
  fft = 7.226362*pow(((15.34432*pow((hour/6.0),1.24))/(51.298686+15.34432*pow((hour/6.0),1.24))),2.520435);
  */
  
  long imat;
  double q,e_0,e_t,fft;
  
  imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
  //creep coefficient
  q = Mm->creep_matstiffchange (ipp,im,ido);
  //asympotic elastic modulus 
  e_0=Mm->eliso[imat].e;
  e_t = q*e_0;
  
  fft = e_t*ft_ratio;
  
  return (fft);
} 


/** 
   function returns actual compression strenght for ordinary concrete
   by Stemberk P. and Tsubaki T.
  
*/ 
double b3mat::give_actual_fc (long ipp,long im,long ido) 
{ 
  /*  
      double ffc,hour;
      
      hour = (Mp->timecon.endtime () - Mp->time)/3600.0;
      
      ffc = 50.58454*pow(((15.34432*pow((hour/6.0),1.24))/(51.298686+15.34432*pow((hour/6.0),1.24))),2.520435);
  */  

  
  long imat;
  double q,e_0,e_t,ffc;
  
  imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
  //creep coefficient
  q = Mm->creep_matstiffchange (ipp,im,ido);
  //asympotic elastic modulus 
  e_0=Mm->eliso[imat].e;
  e_t = q*e_0;
  
  ffc = e_t/1000.0;


  return (ffc);
}