creep_rspec.cpp

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

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

  if(Mm->tempr != NULL)
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+1)] = Mm->tempr[ipp];
  else
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+1)] = 0.0;
}


/**
   function returns time of concrete (end of concrete casting) in days

   @param ipp - number of integration point

   TKr, 3.6.2005
*/
double rspecmat::give_tb_time(long ipp)
{
  double tb;
  
  tb = tb_time/86400.0;//days

  return(tb);
}



/**
   function returns time when drying begins in days

   @param ipp - number of integration point

   TKr, 3.6.2005
*/
double rspecmat::give_th_time(long ipp)
{
  double th;
  
  th = th_time/86400.0;//days

  return(th);
}



/**
   function returns number of times in which complinace function is approximated

   @param ipp - number of integration point

   TKr, 3.6.2005
*/
long rspecmat::give_napproxtime(long ipp)
{
  return(napproxtime);
}


void rspecmat::nlstresses (long ipp,long im,long ido)
{
  long i,nm, ncompstr = Mm->ip[ipp].ncompstr, ncompo;
  vector sigback(ncompstr);

  if (Mp->phase==1)
  {
    // computation of increments of stresses due to the irreversible strain increment
    Mm->creep_phase1(ipp,im,ido);

    // computation of increments of stresses due to the temperature strain increment
    if ((im == 0) && (Mm->ip[ipp].hmt & 1))
    {
      nm=Mm->ip[ipp].nm-1;
      ncompo = Mm->givencompeqother(ipp, 0);
      ncompo -= Mm->givencompeqother(ipp, nm);
      for (i=0; i < ncompstr; i++)
        sigback[i] = Mm->ip[ipp].stress[i];
      Mm->computenlstresses(ipp, nm, ncompo);
      for (i=0; i < ncompstr; i++)
        Mm->ip[ipp].stress[i] += sigback[i];
    }
  }
  
  if (Mp->phase==2)
  {
    Mm->creep_phase2(ipp,im,ido);

    if ((im == 0) && (Mm->ip[ipp].hmt & 1))
    {
      nm=Mm->ip[ipp].nm-1;
      ncompo = Mm->givencompeqother(ipp, 0);
      ncompo -= Mm->givencompeqother(ipp, nm);
      // actualization of previous temperature
      Mm->computenlstresses(ipp, nm, ncompo);
    }
  }
}




/**
   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 rspecmat::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_h_old,deps_a,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_dt - 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;//zatim??!!

    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;

    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 rspecmat::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 into other

  @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 rspecmat::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 from other

  @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 rspecmat::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 ip
   @param compother - number of other components
   @param ipp - first integration point on element
*/

double rspecmat::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 rspecmat::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 rspecmat::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 rspecmat::give_actual_ft (long ipp,long im,long ido) 
{ 
  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);
} 


double rspecmat::give_actual_fc (long ipp,long im,long ido) 
{ 
  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);
}