creep_b3.cpp

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

/**
   function stores increments of total stresses into other 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 b3mat::store_dstresses_add(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];
    Mm->ip[ipp].eqother[ido+(nc)+i] = Mm->ip[ipp].eqother[ido+(nc)+i] + dsigma[i];
  }
}




/**
   function returns total strains from other

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

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

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



/**
   function stores total strains into other

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

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

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


/**
   function returns increments of total strains from other

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

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

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



/**
   function stores increments of total strains into other

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

   TKr, 3.6.2005
*/
void b3mat::store_dstrains(vector deps,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)+i] = deps[i];
  }
}

 
/**
   function returns increments of irreversible creep strains into other

   @param deps_cr - vector of increments of strains
   @param nc      - number of components
   @param ipp     - number of integration point

   TKr, 3.6.2005
*/
void b3mat::give_creepdstrains(vector &deps_cr,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<nc;i++){
    deps_cr[i] = Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times+nc*n_ret_times)+i];
  }
}


/**
   function stores increments of irreversible creep strains into other

   @param deps_cr - vector of increments of strains
   @param nc      - number of components
   @param ipp     - number of integration point

   TKr, 3.6.2005
*/
void b3mat::store_creepdstrains(vector deps_cr,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<nc;i++){
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times+nc*n_ret_times)+i] = deps_cr[i];
  }
}


/**
   function returns increments of irreversible shrinkage strains into other

   @param deps_sh - vector of increments of strains
   @param nc      - number of components
   @param ipp     - number of integration point

   TKr, 3.6.2005
*/
void b3mat::give_irrdstrains(vector &deps_sh,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<nc;i++){
    deps_sh[i] = Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times+nc*n_ret_times+nc)+i];
  }
}


/**
   function stores increments of irreversible shrinkage strains into other

   @param deps_sh - vector of increments of strains
   @param nc      - number of components
   @param ipp     - number of integration point

   TKr, 3.6.2005
*/
void b3mat::store_irrdstrains(vector deps_sh,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<nc;i++){
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times+nc*n_ret_times+nc)+i] = deps_sh[i];
  }
}

/**
   function returns increments of irreversible stress-induced shrinkage strains into other

   @param deps_ss - vector of increments of strains
   @param nc      - number of components
   @param ipp     - number of integration point

   TKr, 3.6.2005
*/
void b3mat::give_stressirrdstrains(vector &deps_ss,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<nc;i++){
    deps_ss[i] = Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times+nc*n_ret_times+nc+nc)+i];
  }
}


/**
   function stores increments of irreversible stress-induced shrinkage strains into other

   @param deps_ss - vector of increments of strains
   @param nc      - number of components
   @param ipp     - number of integration point

   TKr, 3.6.2005
*/
void b3mat::store_stressirrdstrains(vector deps_ss,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<nc;i++){
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times+nc*n_ret_times+nc+nc)+i] = deps_ss[i];
  }
}


/**
   function returns free shrinkage

   @param eps_sh  - shrinkage strains
   @param nc      - number of components
   @param ipp     - number of integration point

   TKr, 19.6.2007
*/
double b3mat::give_shrinkage(long ipp,long ido)
{
  long i;
  long n_ret_times;
  long nc;
  double eps_sh;

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

  return eps_sh;
}


/**
   function stores free shrinkage

   @param eps_sh  - shrinkage strains
   @param nc      - number of components
   @param ipp     - number of integration point

   TKr, 19.6.2007
*/
void b3mat::store_shrinkage(double eps_sh,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<nc;i++){
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc+2+2+n_ret_times+nc*n_ret_times+nc+nc+nc+nc)] = eps_sh;
  }
}


/**
   function stores actual humidity

   @param ipp     - number of integration point

   TKr, 3.6.2005
*/
void b3mat::store_hum(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];
  else
    Mm->ip[ipp].eqother[ido+(nc+nc+nc+nc)] = 0.0;
}



/**
   function stores actual temperature

   @param ipp     - number of integration point

   TKr, 3.6.2005
*/
void b3mat::store_temp(long ipp,long ido)
{
  long nc;

  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 b3mat::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 b3mat::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 b3mat::give_napproxtime(long ipp)
{
  return(napproxtime);
}


void b3mat::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 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, 3.6.2005
*/
double b3mat::b3_law (double t0, double tl, double t, long ipp,long ido)
{
  double n,c0,cd,ac,ag,m,z,r,qf,q,eps_shinf,tau,ht,htl,st,stl;
  double hum,temp,hum_prev,temp_prev;
  double jt,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 b3_law (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 b3_law (file %s, line %d).\n",__FILE__,__LINE__);
    abort();
  }

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

  if (type_e == 1)
    q1=600000.0/e28;//measured
  else